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NUMERICAL ANALYSIS AND FORTRAN PROGRAM FOR THE COMPUTATION 
OF THE TURBULENT WAKES OF TURBOMACHINERY ROTOR BLADES, 
ISOLATED AIRFOILS AND CASCADE OF AIRFOILS 

C. Hah and B. Lakshminarayana 

Department of Aerospace Engineering 
The Pennsylvania State University 
University Park, Pennsylvania 16802 


SUMMARY 

Turbulent wakes of turbomachinery rotor blades, isolated airfoils, and 
a cascade of airfoils were Investigated both numerically and experimentally. 

The study was confined to low subsonic and incompressible wake flows. A finite 
difference procedure was employed in the numerical analysis utilizing the 
continuity, momentum, and turbulence closure equations in the rotating, curvi- 
linear, and non-orthogonal coordinate system. A non-orthogonal curvilinear 
coordinate system was developed to improve the accuracy and efficiency of the 
numerical calculation. 

Three turbxilence models were employed to obtain closure of the governing 
equations. The first model was comprised of transport equations for the 
turbulent kinetic energy and the rate of energy dissipation, and the second 
and third models were comprised of equations for the rate of turbulent kinetic 
energy dissipation and Reynolds stresses, respectively. The second model 
handles the convection and diffusion terms in the Reynolds stress transport 
equation collectively, while the third model handles them Individually. All 
three models were modified for the effect of streamline curvature. The rotation- 
originated redistribution terms were added in the transport equation of 
Reynolds stresses of the second and third models. This is the first reported 
attempt to modify Reynolds stress model for the effect of the streamline curva- 
ture and rotation. The turbulent wakes of an isolated airfoil and a cascade 
of airfoils are handled as simpler cases of the general rotating three- 
dimensional wake. The numerical results demonstrate .that the second and third 
models provide accurate predictions, but the computer time and memory storage 
can be considerably saved with the second model. The experimental data are 
utilized to compare various turbulence closure models for the effect of 
streamline curvature. 
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INTRODUCTION 

Tn the study of modem aerodynamics, turbulence behind turbomachinery 
rotor blades is a significant factor in the establishment of Improved design 
criteria. Accurate information on the characteristics of this flow field is 
also Indispensable in predicting noise level and determining vibration 
characteristics . 

Several experimental investigations and approximate analyses have contrib- 
uted to the understanding of these flow phenomena, although there has been no 
theoretical analysis. The rotor wake is highly complex and three-dimensional 
due to centrifugal and Coriolis forces in the shear layers; blade geometry, 
density, and pressure gradients in radial direction; and other effects. 

Velocity and pressure fluctuate periodically as well as randomly. In addition, 
the turbulence structure is affected by streamline curvature and machine 
rotation. 

The rotor wake has different characteristics as compared to the wake of 
a single airfoil or a cascade of airfoils. Its velocity profile and decay 
characteristics vary in all the three coordinate directions. The analytical 
technique and correlation for a single airfoil cannot be used for prediction 
of the rotor wake, and the cascade does not represent the rotor wake's actual 
flow situation. 

Direct measurement of the rotor wake could provide the data necessary for 
a theoretical analysis of the flow field, but is very expensive and time- 
consuming, Also, there are many parameters which influence the development 
of the rotor wake; blade geometry, solidity, hub-to-tip ratio, speed of 
rotation, inlet turbulence level, three-dimensional pressure gradient, and 
inlet swirl. Therefore, the correlation method of modelling the wake can be 
used only if ample experimental data, including all the important parameters, 
are available. 

The numerical scheme with a suitable turbulent flow model can be con- 
sidered as another promising way of analyzing the rotor wake. However, the 
literature survey (described later) on numerical analysis of two- and three- 
dimensional turbulent wakes indicates that most of the efforts have been con- 
centrated on simple, two-dimensional flow fields without streamline curvature 
and rotation. 

To provide the turbomachinery aerodynamicist or acoustician with an 
accurate method of determining the rotor wake development, or in understanding 
the complex features of the flow, a fully numerical scheme which utilizes 
suitable turbulence closure models is developed in this report. Also, its 
correlation with all the existing rotor wake data and related approximate 
analysis is presented. 


Objectives and Method of Investigation 

The following objectives were pursued during the investigation of the 
turbulent rotor wake reported in this paper: 
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1. To develop 3^ numerical scheme that can predict the characteristics of 
the rotor wake accurately, with reasonable computing time and computer 
memory storage. 

2. To establish a proper turbulence closure model that accounts for the 
effects of streamline curvature and rotation. 

3. To predict the development of turbulent wakes of a single airfoil and 
cascade of airfoils with the developed numerical scheme. 

The problem of predicting the flow field behind rotor blades is formulated 
using the equations of continuity, mean momentum, and turbulence closure. 

In the analysis of complex flow problems, the judicious choice of the 
coordinate system and simplification in the governing equations are pre- 
requisite for successful results. For the present problem, a curvilinear 
coordinate system, which includes strearawise, normal, and radial directions, 
is employed. Though this coordinate system is non-orthogonal, it provides 
simplicity in governing equations and in describing the boundary condition. 

The existing turbulence closure models are critically reviewed for their 
performance in the presence of streamline curvature and rotation. Three 
typical turbulence closure models are considered for the present analysis . 

The first model is comprised of transport equations for the turbulent kinetic 
energy and the rate of energy dissipation. The second model makes use of the 
equation for the rate of turbulent kinetic energy dissipation and the Reynolds 
stress equation; however, the effects of convection and diffusion in the 
Reynolds stress transport equation were handled collectively. The third model 
utilizes equations of turbulent kinetic energy dissipation and Reynolds 
stresses in nearly exact form. All three models are modified for the effect 
of the streamline curvature. The effect of rotation is Included in the second 
and third models through the rotation-originated redistribution term in the 
transport equation of Reynolds stress. 

Bruce Reynolds and A. Ravlndranath made available their experimental data 
for comparison purposes. 
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PREVIOUS INVESTIGATIONS 


The term ”wake*' Is commonly used to refer to the flow region behind a 
solid body, where the effect of viscosity cannot be neglected. Earlier 
investigations were concentrated on the wake at a low Re 3 molds number (laminar 
wake) of simple bodies (flat-plate, cylinder, and single airfoil). With the 
development of turbulence theories in fluid mechanics, the wake at a high 
Reynolds number (turbulent wake) has been widely studied for scientific and/or 
engineering purposes. One of the main differences between laminar and turbulent 
wakes is that the turbulent wake cannot be analyzed with the property of a fluid 
(kinematic viscosity). Hence, flow (turbulence closure) modelling is necessary 
for the analysis of the turbulent wake. To analyze the turbulent near wake, 
the scientist prefers using inner scales based on a turbulence structure at the 
near wake region, while the practical numerical analyst prefers an outer scale 
based on the far wake. This difference in approach will be discussed further 
later in this chapter. 

In the following sections, previous work on analytical approaches, 
turbulence closure modelling, numerical approaches, and experimental investiga- 
tion of two- and three-dimensional turbulent wakes will be comprehensively 
reviewed. Also a review of studies on the wakes of a single airfoil and a 
rotor blade by Raj and Lakshminarayana [1] and Re 3 molds and Lakshminarayana [2] 
and on the wake of a cascade inlet guide vane, and stator wake by Ravindra- 
nath and Lakshminarayana [3] will be given. 


Analytical Approaches 

The turbulent flow field behind a solid body can be divided into very near 
wake, near wake, and far wake regions depending on the evaluation of the flow. 
The far wake region, where the flow evolution is very slow, has been widely 
investigated experimentally and theoretically. The similarity solution exists 
at this wake region because of the slowness in evolution. Since it can be 
found in any standard text (Tennekes and Lumley [4]), the analytical procedure 
for the two-dimensional turbulent far wake is not reviewed here. 

Analytical solutions for the wake downstream of the trailing edge of the 
flat plate and the isolated airfoil have been obtained for a laminar initial 
boundary layer with laminar viscous stresses in the near wake (Goldstein [5], 
Stewartson [6]). For the laminar near wake, the coordinate expansion solution 
for the inner and outer layers can be easily obtained because the flow structure 
is relatively simple compared to a turbulent wake. 

Though numerical handling of the very near wake exists for the simple two- 
dimensional turbulent wake without streamline curvature, no complete analytical 
solution for this two-dimensional wake has been obtained. Robinson [7] 
obtained an analytical solution for the near wake region of a flat plate with 
various assumptions. His solution for the wake centerline mean-velocity defect 
is as follows: 
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where g(Ang - 1) ■ K^Cr^)* is the friction velocity,, and k, Yr» Bj., are 
constants. Although this solution is based on various assxunptions on the flow 
structure at the near wake, it can be further refined and expanded for the 
general turbulent wake. 

The analytical solution for the far wake of a rotor blade was derived 
by Lakshminarayana [8] and Raj and Lakshminarayana [9]. Lakshminarayana [8] 
used a cylindrical coordinate system and provided the wake centerline velocity 
defect and the wake width growth. Raj and Lakshminarayana [9] adopted a 
curvilinear coordinate system for their analytical solution of rotor wake for 
downstream. They Indicated qualitatively the dependency of rotor wake charac- 
teristics on rotor geometry, turbulence quantity, operating condition of the 
rotor, and the streamline curvature. The flow characteristics of the very near 
wake and the near wake of a rotor blade are not well understood because the 
flow evolution is very fast and multiple-length scales, including the stream- 
wise distance from the blade, strongly affect the flow development. The 
governing equations cannot be simplified, and remain non-linear in these 
flow regions. 

Due to the complexity of the flow structure at the very near and the near 
wake of the rotor blade, no analytical or numerical solution has been tried 
before this study. 


Turbulence Closure Modelling 

With the advent of high speed digital computers, the prediction of a 
complex flow field, which cannot be analytically solved, has become a problem 
of major concern. Most interesting flow fields occur at high Reynolds 
numbers and contain random fluctuating components of flow characteristics. 

This results in more unknowns than would be the case for equations in the 
customary description of turbulence. To close the gap between the unknowns 
and the equations, various turbulence closure models have been studied. 

The zero equation closure model, which uses no additional partial differen- 
tial equation for the turbulence closure, has been widely studied and still 
applied for sophisticated practical problems. This model relates the Reynolds 
shear stress directly to the mean shear rate of the flow and cannot be used 
for the flow where multiple length and velocity dominate the problem. Prandtl 
[10] provided the earliest form of this model. Cebecl and Smith [11] have 
suggested modified forms to handle more complex turbulent flows . Their model 
is primarily for two-dimensional flow where only the shear stress component of 
the Reynolds stress tensor is Important. 

In an effort to reduce the empiricism required in closure modelling, 
models Incorporating one or two partial differential equations for the transport 
of turbulence structural characteristics have been developed. Norris and 
Reynolds [12] proposed a model with the turbulent kinetic energy equation, and 
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Wolfs tein [13] introduced a model with the transport equation for the turbulence 
length scale. These one-equation closure models describe the turbulent flow 
field more accurately than zero-equation model, but the application is 
restricted for the two-dimensional flow. 

A two-equation model, which incorporates additional equations for the 
velocity scale and scalar dissipation, has been successfully applied for simple 
two-dimensional turbulent flows. Jones and Launder [14], Launder and 
Spalding [15], and Saffmann and Wilcox [16] have proposed their own versions 
of the two-equation model. Compared to the zero- and one-equation models, the 
two-equation model has less empiricism and can be used for a wider range of 
turbulent flows. Since the additional equation for scalar energy dissipation 
is derived on the basis of isotropic dissipation, the two-equation closure 
model does not predict Individual intensity components of the Reynolds stress 
tensor. Even though the two-equation closure model was successful in predicting 
some simple three-dimensional turbulent flows, general non- isotropic three- 
dimensional turbulent flow fields cannot be described with two turbulence 
structural quantities. Mellor and Herring [17] provided an overview of one- 
and two-equation models, and Reynolds [18] gave a more detailed review of the 
two-equation model as well as other models. No major modification of the two- 
equation model for the effects of streamline curvature and rotation has been 
reported before this report. 

The Reynolds stress model, which utilizes the transport equations of all 
the Reynolds stress components, has been investigated by Hanjalic and 
Launder [19], Launder, Reece, and Rodi [20], and Lumley and Khajeh-Nouri [21]. 
This model has been applied for relatively simple, thin shear layers. The 
numerical predictions using these models have been very promising. Because 
all the non-zero Reynolds stress components are individually estimated through 
the transport equations, this model is more suitable for examining the complex 
three-dimensional turbulent flow. For the computation of turbulent flow fields, 
where multiple-velocity and length scales dominate the turbulent structure, the 
accurate estimation of each Reynolds stress component is essential. In this 
situation, the Re 3 molds stress model works better than the one- or two-equation 
closure models. 

Large eddy simulation, which numerically calculates the three-dimensional, 
time-dependent, large-eddy structure of the turbulence in addition to conven- 
tional modelling of the small-scale eddy, has been under intense development 
by the group at Stanford University. Large-eddy simulations are still in their 
infancy and, therefore, can hardly be used for the calculation of realistic 
problems because of the huge computing time and memory capacity required. 

For the present investigation, the Reynolds stress model was modified for 
the effects of streamline curvature and rotation. This is the first reported 
attempt to modify the Reynolds stress model in this application. Also, the 
two-equation model was applied, and the numerical results are presented for 
comparison purposes . 
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Turbulence Closure Modelling for the Effect of Streamline Curvature 

The effect of streamline curvature in a turbulent flow was first recognized 
by Prandtl [10], and many investigations have been carried out since then. 
Experimental studies of curved botmdary layers have revealed that the turbulence 
structure as well as the mean velocity profile is changed considerably due to 
the effect of the streamline curvature. 

Prandtl [10] demonstrated the effect of streamline curvature with the 
motion of a disturbed element of fluid in an inviscid and Incompressible flow. 
According to his theory, this curvature stabilizes the flow field and, con- 
sequently, diminishes the turbulent transportation, where the angular momentum 
increases with the radius of streamline curvature and vice versa. 

Quantitatively, the following expression can be derived for the criterion 
of stability due to streamline curvature (Bradshaw, [22]): 

where U = tangential velocity, r = radius fo curvature, and w = corresponding 
frequency. A real value of U) means stable, and an Imaginary value means 
unstable. 

Various proposals to Include the effect of streamline curvature in the 
modelling of turbulent flow have been made. Prandtl [10] suggested a modifica- 
tion of mixing-length theory based on his argument of the qualitative effect 
of streamline curvature, A velocity scale was proposed by So [23] for curved 
shear flow from the simplification of Re 3 molds stress transport equations on 
a curvilinear coordinate system. Bradshaw [22] also Incorporated the effect 
of streamline curvature by using the mixing-length concept and the analogy 
between buoyancy and curvature effects, in addition. Launder, Pridden, and 
Sharma [24] modified the two-equation closure model for this factor, and 
Chambers and Wilcox [25] added a singular production term in the equation of 
the turbulent kinetic energy to represent streamline curvature effects. 

All of these approaches are based on the experimental phenomena and bear 
the effect of streamline curvature in various closure models qualitatively. 

The zero-equation model and any curvature modification based on the zero- 
equation model cannot be applied for the calculation of the rotor wake because 
the scaling laws of the turbulent rotor wake are not previously known. Also, 
the application of a simple modification in the two-equation turbulence closure 
model to the three-dimensional rotor wake is in doubt because the turbulence 
structure of the rotor wake is non-isotroplc, 

A more rational approach for the modelling of curvature effect in the 
turbulence closure scheme is pursued in this thesis, and the results are 
compared systematically. 
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Turbtilence Closure Modelling for the Effect of Rotation 


Research literature discussing the effects of rotation on flow fields is 
quite extensive, but very few theoretical investigations have been made that 
fully account for this influence. It was first recognized by Taylor [26] 
during his study of the stability of rotating flow. Several years later, 
Trefethen [27], during his investigation of the flow in a rotating long tube, 
recognized that the transverse pressure gradients by Coriolis acceleration 
could affect the process of transition by a similar mechanism of density 
stratification in a plane flow. 


Bradshaw [28] proposed a local stability parameter for plane-rotating 
shear layers as a gradient Richardson number. 


i rdU^^ 

'^dy- 


es) 


where = rotational speed and R^ > 0 indicates the tendency for rotation to 
stabilize the flow. The above parameter was derived from the analogy between 
flows with mean streamline curvature and horizontally stratified shear layers 
with vertical density gradients, where centrifugal or buoyancy forces produce 
normal pressure gradients. A comprehensive review of the effect of rotation 
on the turbulent boundary layer was given by Johnston [29]. 


Most numerical studies (Launder, Pridden, and Sharma [24]) did not consider 
the effect of rotation on the turbulence structure in the computation of rotat- 
ing turbulent flow. Although the total kinetic energy of turbulence is not 
changed by rotation, the energy is redistributed by rotation when the turbulence 
is not isotropic. Majumdar, Pratrap, and Spalding [30] used the conventional 
two-equation closure model for the calculation of the flow in a rotating duct 
and reported that their predictions were not quite satisfactory at high rota- 
tional speed. 


The effect of rotation is substantial in the turbulent wake of rotor 
blades because the rotational speed is high and the turbulence structure is 
highly non-isotropic. Experimental investigations by Raj and Lakshminarayana [1], 
Reynolds et al [31], and Ravindranath and Lakshminarayana [32] have shown very 
high radial intensity in the rotor wake. This high radial intensity is due to 
the effect of rotation on the turbulence structure, Raj and Lumley [33] 
derived the rlative magnitudes of the turbulence quantities of the rotor wake. 
However, they examined a rotor wake whose streamwise direction is almost parallel 
to the machine axis, and their results [u^ > do not agree with other 

experimental results for general configurations of the wake. A qualitative 
analysis of the effect of rotation on the development of the rotor wake was 
also given by Lakshminarayana and Reynolds [24]. 
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Numferlcal Approaches 

Computational fluid mechanics has rapidly developed into an important 
branch of fluid mechanics with the advent of faster and large computers. 
Numerical approaches have been used primarily for the prediction of complex 
flow fields where other methods do not provide reliable predictions. 

To keep the numerical scheme within the capability of present-day com- 
puters, various numerical schemes have bene suggested. Patankar and Spalding 
[35] proposed a parabolic marching technique for three-dimensional flow. 

Their aim was to solve the three-dimensional flow field in a two-dimensional 
memory array and get the solution through a single sweep in the streamwise 
direction. This technique has been successful for the flow where the stream- 
wise pressure gradient is small and no recirculating zone is involved. However, 
subsequent users have pointed out that an accurate pressure field cannot be 
obtained with the parabolic marching technique in both internal and external 
flows. For the complicated flow field (rotor wake, end wall flow in the 
compressor, etc.), the basic assumption of this technique, that the Reynolds 
stress acting on the plane normal to the streamwise direction is negligibly 
small, does not prevail. Ghia and Studerus [36] and Brilley [37] also presented 
their modified version of the parabolic marching technique; the basic feature 
of the technique is the same and the limitation of the technique is not removed. 

The boundary layer assumption cannot be applied for the rotor wake 
because the streamwise velocity gradient is of the same order of normal velocity 
gradient in the near wake. Therefore, the governing equation should be solved 
in elliptic form to account for the streamwise diffusion of the flow quantities 
at the near wake region. 

The governing equations can be numerically represented with the finite 
element or finite difference scheme. Even though the finite element scheme has 
some merit over the finite difference scheme for the flow whose boundary has 
complex geometry (Zienkiewicz [38]), its application for complex turbulent flow 
problems has not been extensive. Since the boundary geometry and boundary con- 
dition for the present rotor wake problem can be easily described with the 
employed curvilinear coordinate system, the implicit finite difference scheme 
was used. 

To increase the efficiency and accuracy of the finite difference scheme for 
complex flow problems, the body-fitted coordinate system has been widely used 
(Thames, Thomson, and Mastin [39], Sorenson and Steger [40], and others). 
However, this coordinate system has been mostly applied for two-dimensional 
flow problems because huge computing time and memory storage is needed even 
for the generation of coordinate system itself (for example, about an hour 
for simple 3-D geometry on CDC-7600). The non-orthogonal coordinate system 
adopted for the present investigation provides the most merits of the body- 
fitted coordinate system. Also, physically important curvature terms are 
explicitly shown in the curvilinear coordinate system. 

The finite difference equation can be solved by the alternating direction 
implicit (ADI) or successive over relaxation (SOR) methods. The ADI method for 
elliptic problems was first proposed by Peacemin and Rachford [41] and has been 
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widely applied for time-marching solutions for Navier-Stokes equations. For 
the present problem* the SOR method was used because the programming was 
relatively simple and the time-marching solution was not applied for the 
subsonic rotor wake. 


Experimental Investigations 

Properly correlated experimental data would provide a relatively accurate 
estimation of the rotor blade flow field. In previous research endeavors, 
many attempts have been made to find correlations for these turbulent wakes. 
Spence [42] derived the wake centerline velocity decay as follows, from 
systematic measurements of turbulent wakes of a single airfoil: 

1 - ^ = 0.1265 (I + 0.025)"^^^ (4) 

Uq c 

where C = chord at the airfoil, = wake center velocity, = free stream 
velocity, and Z = the strearawise distance from the trailing edge. 

From their experimental data of the rotor wake, Raj and Lakshminarayana [9] 
confirmed that the streainwise velocity component has a Gaussian similarity 
profile at the far wake when proper length and velocity scales are Introduced. 

In addition, they provided a correlation for the decay of streamwise and radial 
velocity components. Reynolds and Lakshminarayana [2] also confirmed with their 
measurements of fan rotor wakes that the similarity profile exists in axial 
and tangential velocity components. However, he found that the similarity in 
the velocity profile is not evident at the very near wake region. He correlated 
all the existing turbulent wake data with the following expression: 


^s - ^c _ ^1/4 



(5) 


where = wake center velocity, Uqo = free stream velocity, S^/C = virtual 
origin of wake, C = chord of blade, S = streamwise distance from the trailing 
edge of the blade, and and are constants which have different values 

for the single airfoil, cascade, and rotor blade wakes. The experimental 
surveys of rotor wake at The Pennsylvania State University (Raj and Lakshmi- 
narayana [9], Reynolds et al. [31], and Ravindranath and Lakshminarayana [32]) 
include the measurements at the near wake region. More parameters are system- 
atically accounted for in their investigations than in other experimental 
studies (Evans [43], Schmidt and Oklishi [44], Gallus [45], and Kool et al [46]. 
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NUMERICAL ANALYSIS OF TURBULENT WAKES OF AXIAL 
FLOW TURBOMACmNERY ROTOR BLADES 


Physical Nature of the Rotor Wake 

Unlike single airfoil or cascade wakes, the rotor wake Is three- 
dimensional (see Figure 1). The Inbalance of the radial pressure gradient 
and the centrifugal forces In the radial direction inside the shear layers 
result in the radial velocity component. Furthermore, the geometry of the 
rotor blade and the density gradients in radial direction also contribute In 
generating these components. Due to the nature of these forces, therefore, 
the radial velocity component has an outward or Inward direction. 

As can be seen In Figure 1, the rotor wake develops In approximately a 
helix shape on the cylindrical surface (s direction) , with 3 as the angle 
between the machine axis and the streamwlse direction. Because a rotor blade 
Is twisted In many cases, the outlet angle from the rotor blade and angle 3 
both change In radial direction. One major curvature component In the three- 
dimensional rotor wake Is due to this change of 3. 

The rotor wake Is periodically unsteady when viewed from a stationary 
frame of reference, and the absolute tangential velocity components look like 
a jet (the maximum absolute tangential velocity occurs at the blade trailing 
edge). On the other hand, the rotor wake Is steady If viewed from a frame 
rotating with the rotor blades. The rotation vector of the rotor blade has 
a direction parallel to the machine axis, while the streamwlse component of the 
rotation vector varies In the radial direction. 

Due to the loading on the blade, the boundary layer on the suction 
surface Is thicker than that on the pressure surface at the trailing edge of 
the blade, and the rotor wake Is asymmetrical. 

The major components of curvature In the rotor wake is shown in Figure 2. 
rjn is the radius of curvature due to machine radii of the hub and annulus 
walls; rg is the radius of curvature due to the change of 3 in the radial 
direction; and tg is the radius of the streamline curvature on the cylindrical 
surface. The effect of streamline curvature depends on the momentum gradient 
in the direction of the radius of the curvature. For the turbulent wake, the 
fluid particles near the wake centerline are most affected by the curvature 
because the velocity gradient Is large near the wake centerline. The effect 
r^ is substantial when the machine radius is relatively small and the blade is 
substantially twisted in the radial direction, r^ and tg effects are dependent 
on the variation of outlet angles in the radial and streamwlse directions, but 
for most practical rotor blades, the effects of these two curvature components 
are cancelled out because the momentum gradients have opposite signs. 
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Figure 1. Characteristics of 


a Rotor Wake 


Streamline 
direction at 
tip 



Figure 2. Major Components of Curvature In Rotor Wake 
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The Curvilinear Coordlante System for the 
Numerical Analysts of the Rotor Walce 

As mentioned earlier, in the numerical analysis of complex flow problems, 
the judicious choice of the coordinate system and the sinpllfication of the 
governing equations are prerequisites for successful results. The rotor wake 
develops between the annulus and hub walls, which for the present study are 
concentric cylinders. The fluid particles from the trailing edge of the rotor 
blade move approximately in a helical path whose angle varies radially. To 
make the governing equations simple, the streamwise direction (s) is chosen 
as one of the coordinates. Because the hub and annulus walls form the boundary 
surface, the radial direction (r) is also chosen as one of the coordinates. 

The third coordinate direction is chosen normal to the streamwise direction (n) 
on the cylindrical surface (constant radius). This physical coordinate system 
fixed to the rotor blade is shown in Figure 3. 

This coordinate system is non-orthogonal, as will be shown later, but 
provides simplicity in describing governing equations and boundary conditions. 
The angle 3 between the streamwise direction and the machine axis varies 
radially. The variation of 3 in the streamwise and normal directions was not 
considered in the derivation of the coordinate system, since rotor wake 
measurements show very small changes of 3 in both directions. 

To get the metric tensor of this coordinate system, the corresponding 
Cartesian coordinates system is introduced at the axis of the rotor shown in 
Figure 3. The y direction coincides with the axial direction of the rotor. 

Any arbitrary point p can be represented both in Cartesian coordinates (x,y,z) 
and in the curvilinear coordinates (s,n,r). The components of (s,n,r) and 
(x,y,z) are related as follows: 


X 


. /"s sin 3 + n cos 3'^ 
-r Sint 


( 6 ) 


y = s cos 3 “ n sin 3 


(7) 


= r cos 


rs sin 3 + n cos 3 


z 


r 


( 8 ) 
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where s, n, and r are atreanwlse, principal normal, and radial directions, 
respectively. 

The derivatives of the Cartesian coordinates with respect to the curvi- 
linear coordinate system are as follows: 


3s 


_j- o fs sin 6 + n cos 
sin g cos 


(9) 


3x 

3n 


o rs sin B + n cos B>i 

cos p cos ■ — ■ J 


( 10 ) 


3x _ _ sin B -f n cos g’ t 

3r r 


- cos 


I'S sin B + n cos B^ T s sin B + n cos B . z' « . dB" 

I J r ^ + (s cos B - n sin B)-^ 


( 11 ) 


3y _ Q 
- cos B 


( 12 ) 


- sin B 


(13) 


3r 


(-S sin B- n cos B)"— 

^ dr 


(14) 


3z 

3s 


= - sin B sin p - ^ 


(15) 


= - cos 


B sin(- 


s sin B + n cos B 


3z 

3n 


r 


(16) 
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3z _ rs sln3 + n cosB-^ 
_ . cos I ) 


oln (^ sing + n cosBj ^ a slnB + n cosB 


+ (s cosB - n sing}^”] 


(17) 


The fundamental metric tensor is obtained from the following 
relation: 


a 


( C = 1.2,3) 

^ ax an‘^ 


(18) 


where = Cartesian coordinates and » curvilinear coordinates. The 
resulting metric tensor In covariant components (S^j) 


®11 * ^ » ®12 “ ®21 ~ ^ » ®22 * ^ 


* » - s sin B -t- n sinBcosB dg 

®13 ^31~ ~ r * " d7 


_ _ _ s SinBcosB + n cos 6 . dB 

^23 «32“ ; + s d? 


2 o o 2 

s , 2„ 2s 


833 ” 1 + ^ sin'^B SinBcosB “ + sinBcosB 


+ cos^B + cosBsinB - cos^B 

r r 


+ 2 S£ slnS)^ + 

r ^dr Mr-' 


( 19 ) 
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and 


•ij 


1.0 


( 20 ) 


Then the metric tensor In contravariant components (g^^) is obtained 
using the following relation. 


gij ^ Cofactor of in 

l®ij| 


( 21 ) 


The resulting metric tensor in contravariant component (g^^) is, 

2 

* 1 + -^ sin^e + sin^BcosS 
r r 

2 2 

, 2 f cos Bsin B . 2 „ . . rd6^2 \ 

+ n { + - cosBsinB^ + ^ 

^12 _21 2 f sin^gcosB sin^B dB ^ 

® ^ ^ ~ r dr 

+ sn { “ sin^Bcos^B - + n^ r s^ inBcos B cos B 

^2 '^dr^ J r dr-* 

2 

13 ^ 31 ^ s sin B + n sinBcosB . ^ 

* ® * r "dr 

_22 _ , « . .2 f sin^B cos^B 2 . ^ . dB 

g » 1,0 + s { 5 SinBcosB ^ 

* r dr 


+ O' } 
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3 2 

o _ _ rsinB cos cos 3 

® -5 — 


4. *0 

+ — ^ cos 3 

r 


23 32 

g . " g 


ssin3 cos3 + n cos 3 
r 


s 


dr 


33 

g 


1.0 


( 22 ) 


The Christoffel symbols of the first kind and the second kind are 
obtained using following relations. 


r 


ijk 


1 

2 


9 xj 




kl 

g 




(23) 


Then the Christoffel symbols of the first kind are as fellows: 


r *r “F *r »r *r iBn 

Ull 112 122 212 211 131 311 


p *r "“F 

^121 221 ‘222 ‘232 ‘322 ^ 


sin^3 


113 


r ■■ F 
123 213 


sin3cos3 
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^132 “ ^312 “ dr 


9 2fl dB . rdB^Z 

•■aia “ ’’l33 “ ® ■ T dr ® MrJ 


T\ 2 Tl 2 \ 

+ slnBcosB ~ ( 7 ^ ^ ' <ir 


2« 
cos B 


231 * 321 " dr 


ti 2 

r » sinBcosB + ^ 

323 ^2 ^ 


o s 2 s . 2.\ dB 

+ [ ^ cosBsinB - ~ cos B + 7 sin Bj 




d^B 

- - "T 2 ®W 


d B ntMl^ 
^332 ■ - ® Z2 - "WJ 


. I- , 2„ o rtft (sin^B-cos^B) 

2 ^ sinj. + ^ SinBcosB “ + W 


^333 - ® 
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slnBcoeg d g . rdg- \ rd g^ t I 

r ,2 '•drJ ^ . 2^ -> 

dr dr 


*- r 


+ » n f- 3 slnBcosS + |i 

T 


9 9 9 

. , 8lngcosgrdg\2 . (sin g-coa g) d g1 
+ T— ld?J ^ r 

0 0 0 
+ r cos g _ 2slngcoag dg ^ ( cos g-aln g ) i -dg -j 2 

^ 3 r dr r ^dr-* 


2 2 
slngcosg d g . 

r , 2 . 

dr dr 




and the Chrlstoffel symbols of the second kind are 


rii - - -T 

r 


" “ ■% sln^gcosg 


22 


13 




- s 2 2 

— 2 sin geos g 

r 

- sin^g - sin^gcosg ^ 

r"^ r 


8^ r s±n%cosg 2^ 2 „ dg 

_ ^ Bln geos g ^ J 


dr 


*33 

,2 

13 


^ 2s^ ^3^ „ dg 

-r sin g - — j- sin gcosg ^ 

r^ r'^ 


2 2 
- sin^g(l + cos^g) sin^gcosg + 

r"^ 


dr 


( 24 ) 
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^11 * f dr " "T 

^22 * f zf ~ ”T 

2 s . « n d3 8.2* 2- 

r,- - — sinScosg -j T sin $cos 6 

12 r dr 

2 2 

■ ”3 sin^Scofl^B - ^ sinBcosB (1 + cos^B) ^ 
r"* r 

^ sln^BcosB + ^ sin^3(l+2cos^6) ~ + s(^~|} 

y,3 ^ sln^B 

11 " r 


„3 SinBcosB 

r w — ' 

21 r 


22 " r 


sln^B - f SinBcosB || 

^ r 


SinBcosB--^ cos^B ^ 
r 


2 2 
.3 s'^ ._2, , 2 b‘‘ 

*^33 3 sin e + -^ 

r r 


SinBcosB T' 
dr 


(25) 
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For thin shear flow, n/r and Cd3/dr)^ are negligibly small for most rotors. 
Therefore the terms containing n/r and Cd3/dr)^ can be neglected In equa- 
tion (24) and these terms are not shown in equation C25) . 


As was mentioned earlier, this coordinate system is not orthogonal (the 
off-diagonal components in the fundamental metric tensor have non-zero terms). 
However, the boimdary conditions at the hub and annulus walls can be easily 
expressed with this coordinate system. Also, the governing equations can be 
considerably simplified for thin shear layers because the normal velocity 
component becomes negligibly small and streamwise variation of flow quantities 
becomes very small far downstream in the present coordinate system. 

Horlock and Wordsworth [47] proposed an orthogonal curvilinear coordinate 
system for the calculation of the three-dimensional boundary layer on the rotor 
blade, their coordinate system can neither handle the radial variation of angle 
3 and nor be applied for general problems like the rotor wake. Although 
Miyake and Fujita [48] used a non-orthogonal curvilinear coordinate system for 
calculation of the laminar boundary layer of a rotor blade with a specific 
configuration, they did not allow for any curvature component of the blade. 
Consequently, the application of their coordinate system is strictly restricted 
to a very small group of rotor configurations. 


The theoretical analysis of the turbulence structure of the rotor wake 
can be done more easily on the physical coordinate system than on the cylin- 
drical polar coordinate system because this system is well geared to the 
physical nature of the rotor wake and important curvature terms are shown in 
the governing equations. Also, a comparison of the rotor wake flow structure 
with those of single airfoil and cascade wakes can be made directly, because 
the physical coordinate system reduces to the Cartesian coordinates used for 
analysis of these latter two wakes. 


Governing Equations in the Rotating Curvilinear Coordinate System 

The present analysis neglects the effects of Mach number on the turbulent 
flow fields and confines to the incompressible, fully turbulent flow. The 
equations governing the incompressible, fully ‘turbulent flow relative to a 
coordinate system rotating with angular velocity, fij:, is introduced in 
generalized tensor form and in the curvilinear coordinates described in the 
previous section. 

The continuity equation for incompressible flow is 


or 




0 







Z 

u 


= 0 


(26) 


0 


(27) 
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4 4 

where U and u are mean and fluctuating contravarlant velocity components. 
Using Christoffel symbols in equation C25), 


au 


± 



0 



In the present coordinate system. 


(28) 


+ iY + iw . « ^ + iz + ^ 

3s Sn 3r * 3s 3n 3r 


(29) 


Where U, V, and W are streamwlse, normal, and radial contravarlant velocity 
components in the present curvilinear coordinate system and u, v, and w are 
fluctuating contravarlant velocity components. 

The momentum conservation equation for laminar flow is 


j j k 




and for turbulent flow. 


ij 




and 


j i J 3u . 1 k „i 

u u, . ■ u — T + u u r. 

J 3x^ ^ 


3u^J . j k „i 
3xJ "" 


(30) 


(31) 


(32) 


where is the alternating tensor, fij is the angular velocity of the rotor, 

p is the fluid density, and p* is reduced pressure (p* = p - p/2 At 

a high Reynolds number, the last term in equation (31) that represents the 
dissipation due to the molecular viscosity of a fluid can be neglected. 

For thin shear flow in turbomachinery blade wakes. 


V « U , W < U . 


(33) 
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Then the components of equation. (31) in the curvilinear coordinate system 
(Figure 3) are as follows (neglecting lower-order convection terms in the 
mean velocity and all the small curvature terms in Reynolds stresses): 


u -p + V 1^ + w ^11 + ^ rf. + 2UW r?;„ + rj 

9s 9n 3r 11 33 13 22 

+ 2UV ■’* ^23 “ 2J2W*sln6 - - — 

. 12 9P* . 13 9P*, 3 1 a — 3 — .... 

•"s "is:’"® 


3s 


+ V -” + w + u^r?, + w^rL + 2uw r?. + v^r? 


3n 


3r 


11 


33 


13 


22 


+ 2UV r ^2 + 2VW F 23 - 2f2W*cos3 = - 


22 3P* 

® 3n ^ 


23 
23 3P* 


ds 


, 3 — 3 2 3 — 

3r’ - fe'™ 


(35) 


U'|^+V-^ + W“ + 2f2V*cos3 + 2J2U*sin3 + W + UW T:, + V TZ^ 
os 3n 9r 33 13 22 


. OTTT. t. 3 , orrTT t.3 „2 ^3 _ 1, 31 3P* . 32 3P* . 33 3P*, 

+ 2UV r^2 + 2W F23 + u - - -[g -^ + g -g^ + g 


23 


11 


3~ 3— 32 2 _3 — „3 

— uw - — vw - ~ w - u - uv r ^2 


(36) 


where U*, V*, and W* are covariant components and U, V, and W are contra- 
variant components of mean velocity. The solution of these equations gives 
the contravariant components of mean velocity which can be transformed to 
physical ones. 


Turbulence Closure Modelling 
T he Modelling of the Reynolds Stress Transport Equation 

Various turbulence closure models have been proposed to simulate the 
behavior of real turbulent flows. The common feature of a good turbulence model 
is its universality to a wide range of flows and simplicity in its use. Two 
relatively complex turbulence closure models are utilized for the calculation 
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of the rotor wake which develops under the influence of curvature and rotation. 
A two-equation turbulence closure model, which adopts two transport equations 
of turbulence characteristic quantities, has been proposed by Daly and Harlow 
[49] and Lamder and Spalding [15]. Launder, Pridden and Sharma [24] extended 
this model to Include the effects of curvature by modifying the primary decay 
term in the equation for the rate of dissipation. This modification is based 
on the experimental data and approximately represents the effect of curvature. 

To include the effects of curvature and rotation, the Reynolds stress model 
is used in the present analysis. This model, which utilizes approximate trans- 
port equations of Reynolds stress, has been proposed by Rotta [50] and 
Hanjallc and Launder [19]. Raj and Lakshminarayana [1] give the transport 
equation in the rotating curvilinear coordinate system, which is: 


„ i k , k ,,i , 1 . i k . or> k 

U,u u, . + u u.U, . + u u.U, . + u u u., , + 2 Si 2 E u u 
J J J J 1 J J J ^ m 


^kJlm i 

V 



. X 
+ u p 


*k 


) + Vg 


£j 


. k i . i k , 
(u + u u,^ J 


(37) 


where p’ is the fluctuating static pressure and V is the kinematic viscosity. 


To get closure, several authors have suggested a modelled form of the 
above equation. All the models proposed hitherto are for non-rotating 
systems. Most calculations for thin shear flow assume that the convection 
term and diffusion terra in equation (37) are equal. So [23] applied this 
assumption for the calculation of two-dimensional curved flow. For the present 
analysis, the diffusion and convection terms are considered collectively. The 
combined effect of the two terms are assumed to be related to the production 
term. The modelling of the pressure-strain correlation term is based on Rotta 
[50] and Naot, Shavit, and Wolfshtein [51]. Furthermore, the rate of turbulence 
energy production, which appears in the modelling for the mean-strain effects 
in the pressure-strain correlation, is replaced with the rate of dissipation. 
Also, a high local- turbulence Reynolds number is assumed and the scalar represen- 
tation of dissipation (the contraction of the dissipation tensor) is introduced 
for the dissipation term of equation (37). With these models, equation (37) 
reduces to 


0 = (1 + C^) (1 - Y)-2(e=^H^u^uJ + e-^H^u^u^)- |g^^G(l - Y) 

k 


- “(u^uJ - g^^k) 


(38) 


where k = 1/2 g^^u^uJ, £ = the rate of energy dissipation (last term in 
equation 37), and C^x> ^1» "Y constants. 


These equations are still intractable unless some realistic and practical 
assumptions are made. Using the assumptions utilized for the simplification of 
the mean momentum equations, the component form of equation (38) can be 
expressed as follows for the coordinate system of Figure 3. 




0 = 2(1 + C^) (1 - y){[- uv - uw(- sin3cosB + s H"!" + U 


12 
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+ w r 23 > + [-u^ + uw( ^ s±n^3)](|f + u + W rj^)} 

- Afl [-u^ — sln^B + uv(- — sinBcosB + s 4^) + uw (1 
nr r dr 

2 2 2 
+ ^ sin^e " df sinBcosB)] ’ |d + % slnS)e(l - y) 


2 

V : 


■ Si k‘“^ ■ 1 '^ ^ sln'’ 6 )k] 

r 


(39) 


0 = 2(1 + C )(1 - y)[-v^ - vw(- — sinBcosB + s 4^) ] (-^^ + U 

1 r dr dn IZ 

2 . c 0 2 cj Hr 

+ W r_„) + An [uv(- — sin B) + V (- — sinBcosB + s ~) 
s r r dr 

2 2 2 

+ vw(l + sin^B - 2 — 4^ sinBcosB)] - ^{1 + ~ sin^Bcos^B 
^ r dr J l 

r r 

2 — 2 
" dl si"ecosB)c(l - y) “ C ^ sin^Bcos^B 

^ r 

2 

- sinBcosB)] (AO) 


0 = 2(1 + Cj^)(l - y){“VW - w^(- sinBcosB + s ^ 1 (|^ + U 

+ w r^j) + f + « Til + Mr^3)l-4i\[w 

+ (- ~ sinBcosB + s 4^)] + An (uw - — sin^B) 

r dr n r 

- ! - "> - Si kf“' - f > 


(Al) 
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uv: 


0 = (1 + C^)(l - y){[-v^ “ vw(- ^ sinBcose + s |f) ] (|^ + U r ^2 


1 c n p4TT 1 1 

+ w r23> + (-UV + vw ^ sin B) ^ + W 


+ [-UV - uw(- sinBcosB + s 1 (|“ + U 


- 2{fi^[-uv sin^B + v^( - sinBcosB + s 


2 2 — 

- vw(l + sin^B - sinBcosB)] - [-u^ — sin^B 

2 r dr s r 

r 

2 2 

+ uv(- — sinBcosB + s 4^) + uw(l + sin^B - 4^ sinBcosB)]} 

r dr ^ r dr 

r 

2 2 

- ■|(^ sin^BcosB “ ~ f" sin^B)e(l - y) - C —[uv 

r ^ 


2>s^ .3 s^ dB . 2 .. , 

- sin BcosB -r- sin B)k] 

j ^ r dr 

r 


(42) 


vw: 


0 = 


(1 + C^)(l - y){[-v^ - vw(- 4 SinBcosB + s J (|~ + U F^^ 


+ W F 23 ) + (-UV + vw 4 ^ ^^13^ ■*'11“’' 


- w^(- - sinBcosB + s T~) ] (-^ + U F )} - 2{fi [v' 


+ vw(- — sinBcosB + s 4^) ~{ - fi (uv - vw — sin^B) 
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2 

- ^ slnBcosB)]} - sinScosB - s ^)e(l - y) - C 'fl™ 

- einScosB “ s ||)k] (43) 
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0 = (1 + C^)(l - y){[-vw - w^(-^ sinScosB + s 1 (|“ + U 

+ w r 23 > + [-^ 7 ^ ^13 ■*■ 

- 7^(- ^ sin3cos3 + s + u r ^2 + ^ 

+ uw — sin^g) (-^ + U r^_ + W r^_)l - 2\q. [uv + uw(- — singcosg 

r ds ii 13 ‘ s r 

+ s -T^) ] - tu^ - vw(- — singcosB + s 3 ^) - w^(l + 
dr n r dr 

2 2 

+ sin^g - ™ singcosg) ] 1 + "I ^(sin^g)e (1 - y) 

r 

" ^(J)l f ^ r 

where U, V, W, u, v, and w are jaean and fluctuating velocity components in 
s, n, and r (Figure 3) directions. Q and Q are the components of angular 
velocity in s and n directions. ^ ^ 


A Qualitative Discussion of the Effect of Rotation and Curvature 
on the Turbulence Structure 

The effects of rotation on the turbulence structure is qualitatively 
analyzed employing equations (41-44). To enable simplification, the curvature 
terms are neglected. The resulting equations are as follows: 
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where “ (1 + 


( 45 ) 

(46) 

(47) 

(48) 

(49) 

( 50 ) 
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For the coordinate system shown in Figure 3, flg = ficos3, and £2^ =* -fisin3; 
hence the effect of rotation is to decrease i? and increase ^ for negative 
values of uw, while the negative value of w increases the radial component of 
intensity and decreases the normal component of Intensity (^) . These 

effects are reversed when the correlations are positive. The wake measurements 
reported by Lakshminarayana and Reynolds [34] indicate that w is negative. 
Thus, the overall effect of rotation is to increase the radial intensity and 
decrease the streamwise intensity in a compressor rotor wake. Similar 
comments can be made with regard to shear stresses. Negative values of uw and 
W increase the streamwise shear stress C-uv) through the rotation. The radial 
shear stress C“w) Increases when u^ > ^ and when uv is positive. The more 
quantitative estimation of the rotation effect can be made on the basis of 
order analysis. The relative magnitudes of three intensity components are 

u^ > ^ > v^ for a non-rotating flow without curvature. The rotation- related 
terms in equations (45), (46), and (47) contain the term whose 

order is less than 1 for the present problem. The shear stress terms are one 
order smaller than the Intensity terms. Therefore, the rotation-related terms 
in equations (45), (46), and (47) are one order smaller than the intensity 
terms. Also, neglecting smaller coefficients in equations (45), (46), and 
(47), the relative magnitudes of three intensity components can be represented 
approximately as. 


2 (w^) + 4f2 — (-uw - vw) 

w n E 


(u^) + 4fi — uw 

n E 


(51) 


— s (v ) + 4fi — vw 

2 n E 

v = 

(u ) + 4f2 — uw 

n £ 


(52) 


where subscript n represents the value with no rotation effect. 


Equation (51) shows that the relative magnitude of the radial intensity 
as compared to the streamwise intensity is increased by the rotation effect 
for the case of negative values in uw and w. On the other hand, the relative 
magnitude of the normal intensity as compared to the streamwise intensity does 
not change much due to the rotation. Therefore, the relative magnitudes of 

~2 ~~2 ~2 ~2 ”2 

three intensity components become u^w >v or w^u >v , depending 
on the order k/E. Similar estimation can be drawn for the shear stress 
components, but the effect of rotation on them is more complex. Comparison of 
the data for the isolated airfoil (Hah and Lakshminarayana [52]), cascade 
(Raj and Lakshminarayana [53]), and rotor wakes (Lakshminarayana and Reynolds 
[34]) provides a confirmation of the trends predicted by the qualitative 
argument . 
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The major components of the curvature in the rotor wake were shown in 
Figure 2. The curvature components due to the machine radius and the variation 
of the flow outlet angle Cd3/dr) are explicitly included in the governing 
equations of the present curvilinear coordinate system. 

The effect of streamline curvature on the turbulence structure and the 
performance of various turbulence closure models for the curvature effect is 
explained in great detail earlier. 

To represent the effect of streamline curvature more accurately with the 
existing second-order turbulence closure schemes, the transport equation of 
the rate of turbulent kinetic energy dissipation should be modified. For the 
calculation of the rotor wake, two different terms of the energy dissipation 
rate equation were used and the results will be compared. An examination of 
equations (40) through (44) reveals that the effect of the curvature due to 
the machine radius is substantial for a rotor of small radius. The curvature 
effect due to radial variation of the outlet angle depends on 3 and d3/dr. 
Further, this component cannot be neglected for rotors with highly twisted 
blades. Additionally, in the present numerical scheme, the effect of curva- 
ture due to the radius rg in Figure 2 is implicitly included when the modified 
term of the transport equation is used for the rate of turbulent kinetic energy 
dissipation. 


Closure of Turbulence for the Numerical Analysis of the Rotor Wake 


Continuity and momentum conservation equations governing the flow have to 
be solved simultaneously to get the unknown flow quantities (three mean- 
velocity components and static pressure) in the rotor wake. As the mean 
momentum equations include the Reynolds stress terms, turbulence closure 
equations also should be solved to get these stress terms. With the modelled 
transport equations of Re 3 ?nolds stresses, the six components of these stresses 
can be calculated if the values of the turbulent kinetic energy, the energy 
dissipation rate, and the mean shear rates are known at a given location of 
the flow field. 

The transport equations for the turbulent kinetic energy and the rate of 
turbulent kinetic energy dissipation are 


u — - s - c^2 T t ^ 


9x 


e 


k Z ^ 

. r *k 

9x 


(53) 
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(54) 


where P = -Un-u^U,^: C^, , C_o, C_, and C are constants: and S xs the source 
term of the energy dissipation equation. 


The following two forms of the source term in equation (53) were used 
for the cdlculation of the rotor wake: 
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® k “i“j 


( 55 ) 


and 




(56) 


The first form of the source term, equation (55), has been widely used 
for the calculation of thin shear flow (Hanjalic and Launder [19], Launder, 
Reece, and Rodi [20], Launder and Spalding [15], and others), but is considered 
to be mainly responsible for the poor performance in predicting curved 
turbulent flow. The transport equation of the rate of turbulent kinetic energy 
dissipation (which i s d e fin ed as the contraction of the second-order dissipa- 

3ui 3u-i 

tion tensor e,--; = V -;7 — in the Cartesian tensor) for a high Reynolds number, 
« 1 1 oxT, 3x1, 

xs as follows: X- x 




+ U^.(s' u^’^) - v[s'^'"u^.u. . + s'^V.u, . 

J ik *k i,j ’x k,j 


, , i,ki , , i,ik, „ ni£i mkr~; ; 

+ s' U-” u,, + s’ .U-" u, J - 2v g Jg g [sj. s' , 


+ s( s!, - ] - (eu^ ) » . “ “ pi., + si, pin 1 * 

£m,n xk,j *j p ^ £m xk xk^ £m 


(57) 


The above exact form of the turbulent kinetic energy dissipation equation 
does not include the mean strain rate. Furthermore, the use of the production 
of turbulent kinetic energy as the source term in equation (53) does not have 
any theoretical basis and is probably faulty. This point will be further 
discussed later. The second form of the source term, equation (53), is based 
on the term originally suggested by Lumley and Khajeh-Nouri [21], 


S = (uV - I k 6. ,)(u^u'‘‘ - 4 k 6. .) 

3 xj 3 ji 


(58) 


To make the dimension correct, additional terms were multiplied. Since the 
derivation of this terra is based on the order estimation of equation (57) 
and does not include any mean strain rate, it is expected to perform better 
for analyzing the effect of streamline curvature. 

The present turbulence closure model solves two transport equations for 
the turbulence kinetic energy, equation (54), the rate of turbulence dissipa- 
tion, equation (53), and the mean momentum and continuity equations (29, 34-36) 
simultaneously. The Reynolds shear stresses and the components of intensity 
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are calculated through the modelled form of the transport equation of Reynolds 
stress (equations 39-44) . This procedure thus leads to simultaneous solution 
of U, V, W, p, k, and £■ from a coupled set of equations- The components of 

stresses uv, vw, uw) are derived from the Reynolds stress equa- 

tions (uncoupled from the main program, equations 39-44). Though this model 
does not include the individual effects of convection and diffusion, it does 
fully account for the effects of rotation and curvature. Also, the consistency 
of this turbulence closure model is achieved through the coefficient Cj in 
equation (38). Therefore, this model is efficient for the calculation of thin 
shear flow under strong effects of rotation and curvature with considerably 
less confutation time and computer storage compared to the full Reynolds 
stress model. The numerical values for the universal constants used in this 
analysis are given below (Launder et al- [20] and Pope and Whitelaw [54]). 

A slightly higher value of C 01 resulted in a better prediction in the mean 
velocity as well as the turbulence quantity. 


c 

s 

^4)1 

c 

E 



Y 

0.1 

1.5 

0.15 

1.45 

1.9 

0.6 


The value of in equation (56) was 1.8 for all the calculations. 


Numerical Technique 


Solution Scheme 


The continuity and the time-mean momentum equations (29, 34-36), along 
with turbulence model equations, were solved elliptically . At the near wake 
region, the streamwise velocity gradient is of the same order of magnitude 
as the normal velocity gradient; hence, the conventional boundary layer 
calculation cannot predict this region well. According to experimental data 
by Raj and Lakshminarayana [1], Reynolds et al. [31], and Ravindranath and 
Lakshminarayana [32], the Reynolds stresses acting on the nr plane are of 
comparable magnitude to those on the sn or sr plane. This may be due to the 
effects of rotation and curvature. Thus, the parabolic marching technique 
cannot be used. Such strong elliptic effects in the rotating flow was also 
pointed out by Majumdar, Pratrap, and Spalding [30]. The pressure field was 
first assumed to be constant on the plane of constant r and the pressure 
distribution in the radial direction was assumed to be the same as that at the 
free stream. The pressure field was corrected at each iteration using the 
Poisson equation derived from the continuity equation. The procedure is 
similar to those used by Patankar and Spalding [35] and Briley [37], but 
streamwise diffusion was included in the continuity equation. 

The major difference between the present numerical scheme and the so-called 
parabolic marching technique is that the variables downstream are included in 
the finite difference equations and streamwise diffusion of flow quantities are 
included in the numerical solution scheme. 



35 


As will be shown in the following section, the parabolic inarching scheme 
does not predict the wake decay rate accurately compared to the elliptic scheme. 
This is because the parabolic marching scheme neglects the strearawlse diffu- 
sion of flow quantities that is of comparable magnitude to the diffusion in 
the normal direction at the near wake. 


Initial and Boundary Conditions 

As the governing equations are solved in elliptic form, a boundary condi- 
tion on the whole perimeter of domain is necessary. The computing domain 
extended from hub to tip, one complete passage of the wake, and two chords 
downstream of the trailing edge. The domain of calculation is shown in 
Figure 4‘. A, B, C, and D lie on the same cylindrical surface (the annulus 
wall), as do E, F, G, and H (the hub wall). The botjndary curves AD and BC, 
and EH and FG, lie on the same cylindrical surface and make the same angles 
with machine axis. AB and EF are part of helix, and their lengths are the 
same as the spacings of the blades at each corresponding radius. 

Because an elliptic scheme was used to solve the governing equations, 
boundary conditions were specified on all boundary surfaces (ABFE, BCGF, 

DCGH, ADHE, ABCD, and FGHE) . Due to the nature of the rotor wake, a periodic 
boundary condition was applied between the planes ADHE and BCGF (the distribu- 
tion of flow quantities on the surface ADHE is the same as that on the surface 
BDGF) . 

The available experimental data was used for the boundary conditions on 
the surfaces ABFE, DCGH, ABCD, and EFGH, and extrapolation of this data was 
used for the unprovided conditions on the same boundary surface. The boundary 
value of each dependent variable on the boundary surface was obtained from 
the data, except for the rate of turbulent kinetic energy dissipation. 

Because this rate cannot be obtained from the experimental data, it was assumed 
that it Is equal to the production of the turbulent kinetic energy at that 
location. Fine mesh was used in the wake center region and no wall function 
was utilized. 


The Grid System and the Finite-Difference Equations 

The t 3 Tplcal three-dimensional mesh used for the present numerical 
calculation is shown in Figure 5. The turbulent kinetic energy and the rate 
of dissipation were handled as diffusion terms in numerical calculation and 
located at the same point as the streamwise velocity component and static 
pressure. The turbulent Reynolds stress components are computed with the 
turbulent kinetic energy, the rate of dissipation, and the mean velocity 
gradients at given points in the present turbulence closure scheme. Therefore, 
the link between the mean shear rate and the Reynolds stress remains strong. 
Furthermore instability in numerical calculation due to the absence of such 
a strong link (reported by Pope and XThitelaw [54] in their calculation using 
full Reynolds closure scheme) was not sensed in the present calculation of 
rotor wake. The relative dimensions of typical three-dimensional mesh were 
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Figure 4. Calculation Domain 


I 
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As/An = 1.2, Ar/An - 4. Since the radial variation of flow quantities was 
mild compared to the normal and streamwlae variations, relatively large 
distances were used between radial stations. Up to 60, 90, and 20 nodal points 
were used In streamwise, normal, and radial directions. 

As three momentum equations are used for the calculation of three mean- 
velocity components, the continuity equation C29) Is used to calculate the 
static pressure. However, since equation C29) does not Include any pressure 
terms, it should be transformed to a more convenient form. Patankar and 
Spalding [35] and Briley [37] combined the momentum equations for the normal 
and bi-normal velocity components and the continuity equation for the static 
pressure. However, they did not handle the equation in Its exact form and 
simplified the obtained Poisson equation substantially. 

For the present investigation, a Poisson equation for the pressure Is 
obtained by combining the continuity equation and the divergence of the mean 
momenttnn equations. In the rotating curvilinear coordinate system used here, 
the Poisson equation across the incompressible turbulent rotor wake is as 
follows : 

„1 „3 -.2 „3 V „ ijk,« „ V 1 j 

- - U.J e - u 

(59) 

The present rotor wake calculation in this study has neglected the viscous 
term in the equation (59) because the local turbulence Reynolds number is very 
high. The three mean momentum equations (34), (35), and (36) and Poisson 
equation for the pressure distribution (59) were represented in finite difference 
forms. All spaclal derivatives were computed by second-order central differences. 
The velocity component part of the convection term (l.e., of U^U^j) and 
source terms in the momentum equations were evaluated from the values of 
previous iteration for the convenience of the calculation. 

The finite difference forms of the equations (34), (35), (36), and (59) 
are as follows : 
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where Gil = U^-, = |^ + Ur?-, + V r?-- + Wf?-. 
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and generally Gij = U^, = -r^ + U^r^. 
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Sequence of Calctjlatlon, Stability, and Accuracy 

The flow chart of numerical calculation is shown in Figure 6. The 
following steps were used for the calculation. 
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s=s+ds 




Figure 6. Flow Chart of Rotor Wake Code 
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1. The initial distribution of xjnknovm variables were obtained by 
downstrearawiae parabolic marching witb the gJven boundary conditions. 

2. At each new iteration, the Reynolds stress components were confuted 
with the equations C39) through C44) using the previous values of 
mean velocity, turbulent kinetic energy, and the rate of turbulent 
kinetic energy dissipation. 

3. Using the computed Reynolds stress and previous static pressure 
distribution, the momentum equations and two transport equations for 
the turbulent energy and the rate of dissipation were calculated for 
the new values of the three mean-velocity components and the turbulent 
kinetic energy and the rate of dissipation. 

4. New static pressure distribution was obtained using the newly obtained 
mean velocity components and the Reynolds stress components. 

5. The Reynolds stress components were newly evaluated using the new 
values of the mean velocity components, the turbulent kinetic energy, 
and the rate of turbulent kinetic energy dissipation; steps 3 and 4 
were repeated. 

6. The above steps 2-5 were repeated downs treamwise for the whole flow 
field. Since the governing equations are solved in elliptic form, 
the values of dependent variables on the (i-l)th, 1-th, and (l + l)th 
s treamwise stations are used to calculate new values of dependent 
variables at i-th strearawise station. 

7. Whenever the whole flow field is newly Iterated, the convergence 
criteria is tested and if not converged, the whole domain is reiterated. 

The algebraic equations obtained during the numerical steps were solved 
by SOR. Since the present numerical scheme adopts implicit finite-difference 
simulation of the governing equation, no specific source of instability is 
considered. 

Because three-dimensional storage is needed for the elliptic calculation, 
the required memory size was up to 800 K words, which exceeds the available 
computer main core storage. Therefore, an auxiliary memory system was utilized 
for the calculation. 

The computing time per each iteration of the entire field was up to 
four minutes using the IBM 370/3033 at The Pennsylvania State University. 

In Figure 7, the calculated decay rate of wake center strearawise velocity 
is shown at various iterations for the rotor wake data by Ravindranath and 
Lakshmlnarayana [32]. The results in this figure show that the numerical 
scheme gives uniform convergence, and no oscillation is observed. 

For all the calculations reported, the following criterion was used to 
determine the convergence of calculation; 




Figure 7. Wake centerline Velocity Prediction at Different 
Iterations 
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N 





s 0.007 


(64) 


Where U^. = wake center streaniwtse velocity, subscript i = the value at i-th 
strearawlse station, and superscript n = the value at n-th iteration. For 
most calculations, the convergence with the above criterion was achieved at 
the iterations around 20. To check the numerical behavior of computation, the 
streamwise mesh size was halved and the same calculation was done. With the 
smaller mesh size, the convergence was achieved around 14 iterations and no 
overall computing time was reduced. Comparisons between the numerical results 
and experimental data show that the convergence criteria (equation 64) 
provides a sufficiently accurate solution. 
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Table 1. Relevant Data on the Compressor and the Operating Conditions 



Raj and Lak- 
shmlnarayana 
[1] 

Reynolds 
et al. [311 

Ravlndranath 
and Lakshmi- 
narayana [321 

hub /tip 

0.44 

0.44 

0.5 

R.P.M. 

1010 

1010 

1066 

number of blade 

12 

12 

21 

blade shape 

symmetric 

symmetric 

cambered 

ratio of blade chord ( C) 
to spacing (S) at mid- 
radius 

0.68 

0.68 

0.6 

tip radius 

0.27 m 

0.27 m 

0.46 m 

incidence at mid-radius 

0° 

10®, 15® 

7® 

stagger angle at mid-radius 

45® 

45® 

30® 

drag coefficient (Cj)) at 
mid-radius 

0.0055 

0.011 

0.0116 

lift coefficient (C^) at 
mid-radius 

0.028 

0.5 

0.88 

experimental technique 

stationary 
tri-axial 
hot wire 
probe 

stationary 
tri-axial 
hot wire 
probe 

rotating hot* 
wire probe & 
static pres- 
sure tube 

radial stations at which 
measurements was taken 
(r/r^) 

0.488, 

0.536, 

0.581, 

0.6281, 

0.72, 

0.966, 

0.815 

0.86 

0.465, 

0.535, 

0.629, 

0.721, 

0.814, 

0.860, 

0.907 

0.5676, 

0.6581, 

0.7297, 

0.7973, 

0.8615, 

0.9329, 

0.9595 

1 
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COMPARISON BETWEEN ROTOR WAKE DATA AND NUlffiRICAL PREDICTIONS 

Rotor Wake Data and Predictions with 
the Present Turbulence Closure Model 


The present turbulence closure model and numerical scheme were used to 
predict turbomachinery rotor wakes, utilizing the energy dissipation equation 
with the conventional source term, equation C55). The results of the calcula- 
tions using this equation including the modified source term equation (56), 
will be discussed in the following section. 

Three different sets of experimental data of the rotor wake were compared 
with the numerical predictions. Some relevant data on the compressor and the 
operating conditions are given in Table 1. 

For the present investigation, the hub wall and annulus wall regions, 
where no experimental data or flow model were available, were not considered. 
The boundary surfaces ABCD and EFGH in Figure 4 were cylindrical surfaces 
closest to the hub and annulus walls where measurements were taken. The inlet 
boundary surface ABFE and outlet boundary surface DCGH were the surfaces 
closest and farthest, respectively, from the trailing edge of the blade where 
measurements were made. For all three sets of experimental data, the distance 
between the blade trailing edge and the inlet boundary surface was less than 
1/50 chord length of the blade. The outlet boundary surface was located about 
2 chords downstream. 

For elliptic calculation, the inlet boundary condition specified on 
ABFE, and the boundary conditions on ABCD, EFGH, and DCGH, were based on 
experimental data. Due to the nature of the flow field around rotor blades, 
a periodic boundary condition was specified on ADHE and BCGF. Fourth-order 
power laws were used for the representation of measured flow quantities to 
avoid the effect of scatter in the experimental data on the calculation. 


Comparison Between the Numerical Predictions and the Rotor 
Wake Data by Raj and Lakshminarayana [1] 

The comparison between the numerical predictions in this thesis and 
experimental data by Raj and Lakshminarayana [1] is made for the measurements 
at r/r^ = 0.58. Since the experimental data do not include the static pressure 
distribution, the static pressure was assumed to be constant and radial varia- 
tion of the reduced static pressure due to the centrifugal force was considered 
during the calculation. Comparison between the predictions and the experimental 
data for three mean-velocity components are shown in Figure 8. The streamwise 
velocity component U, the normal velocity component V, and the radial velocity 
component W are normalized by Uqo* The turbulence intensity and Reynolds stress 
are plotted in Figures 9, 10, and 11, respectively. In these figures, s is 
the distance from the trailing edge and n is the distance from the wake center, 
n = s = 0 represents the blade trailing edge and S is the blade spacing. The 
prediction is quite accurate for the streamwise component of the mean velocity. 
The discrepancy in the normal velocity component may be due to inaccurate infor- 
mation on the boundary condition or scatter in the experimental data. 
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Figure 8. Profiles of Mean Velocity Components at r/r^-O-SS 
(Data of Raj and Lakshminaray ana , [9] 
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Figure 9. Reynolds Stress at r/rt=0.58 and s/C=0. 335 
(Data o£ Raj and Lakshminarayana , [9] 
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Figure 11. Reynolds Stress at r/rj-=0.58 and s/C=1.069 
(Data of Raj and Lakshminarayana , [9] 




53 


The prediction of the radial velocity component shows slower decay of 
velocity peaks than the experlm^ital data. This may be due to Improper turbu- 
lence closure modelling or experimental error due to the finite distance 
between the trl-axLal probe wires. The prediction of the turbulence intensity 
shows good agreement with the experimental data, while the prediction for the 
Reynolds shear stress shows considerable deviation from the reported results. 
The discrepancy in the Reynolds shear stress also may be due to scatter in the 
experimental data which were obtained with the stationary tri-axial hot-wire 
technique. Howevar, overall agreement between the numerical prediction and 
the experimental data are quite satisfactory if the complexity of the flow 
field is considered. 

The typical distribution of the rate of turbulent kinetic energy dissipa- 
tion across the wake is shown in Figure 12 , 


Comparison Between the Numerical Predictions and the Rotor Wake Data 
by Reynolds et al. [31] 

The second set of data for the comparison were obtained from Reynolds et al. 
[31]. They measured the wake in the same compressor as Raj and Lakshminarayana 
[1], but the rotor was operated with higher blade loadings. Here again, the 
experimental data did not include the static pressure distribution of the whole 
flow field, and this pressure was assumed to be constant. The radial variation 
of the reduced static pressure due to centrifugal force was considered during 
the calculation. 

Comparisons of the mean-velocity, turbulence-intensity, and shear-stress 
profiles at various axial locations and r/r^- = 0.72 are shown in Figures 13 
through 16. Strearawise velocity components were accurately predicted for both 
incidence angles. Even though some minor disagreements existed between the 
prediction and the experimental data, the radial and normal velocity components 
were well predicted. The experimental data showed that the peaks in the radial 
velocity distribution become wider and smaller with the increase of the incidence 
angle, this trend is well predicted by the numerical calculations. For the 
present sets of data (1 = 10®, 15®), the boundary layers at the blade trailing 
edge are thicker than those of the data in the previous section (i = 0®), and 
the effects of rotation and streamline curvature are substantial compared to 
the first example. According to the qualitative analysis on the effect of 
rotation, the radial turbulence intensity is Increased due to the rotation, 
and this is nearly equal to the streamwise intensity for the present experi- 
mental data. Again, this trend is more than adequately or, "successfully 
predicted" by the numerical calculations. Therefore, the present method can 
fully account for the rotation effect through the modelling of the Reynolds 
stress equation in the rotating curvilinear coordinate system. The typical 
distributions of the rate of turbulent kinetic energy distribution are shown 
in Figure 17 for the incidence angles of 10® and 15®. 
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at s/C=0.066 



Figure 12. Energy Dissipation Rate at r/rt=0.58 
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Figure 14(b). Reynolds Stress at r/rt=*0.72,s/C=0.291 and 
i*15® (Data of Reynolds et al., [31]) 
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Figure 15(b). Reynolds Stress at r/rt=0.72, s/C=0.94 

and i=15® (Data of Reynolds et al., [31j) 
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Figure 16(a). Reynolds Stress at r/r^.«0. 72,s/C«0.93 

and i-lcP (Data of Reynolds et al., [31J) 
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Figure 16(b). Reynolds Stress at r/rt=0.72, s/C=1.65 
and 1=15° (Data of Reynolds et al., [31]) 
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Comparison Between the Numerical Prediction and the Rotor Wake Data 
by Ravindranath and Lakshmiiiarayana [32^3] 

In the third comparison of experimental data, the findings of Ravindranath 
and Lakshminarayana [32, 3] were used. They measured the wake in a single-stage 
compressor (Table 1). The rotor was non-zero camber. The experiment was done 
near the design condition, and the blade was moderately loaded. The data 
include a complete survey of the static pressure distribution across the wake 
at mid-radius. For the numerical calculation, the radial variation of the 
reduced static pressure was assumed to be dependent only on the centrifugal 
force. Also, the rate of turbulent kinetic energy dissipation was assumed to 
be equal to the production of turbulent kinetic eneryg. 

The comparison of mean velocity is given in Figure 18. The overall agree- 
ment between the numerical prediction and the experimental data is good, 
although there is some discrepancy in the profiles of radial velocity components. 

The comparisons of Reynolds stresses are given in Figures 19, 20, and 21. 
Even though the numerical results do not agree exactly with the experimental 
data, the overall characteristics in the local distribution of Reynolds stress 
are excellently predicted. The numerical prediction does accurately represent 
the qualitatively analyzed effect of rotation, that the radial turbulence 
intensity is increased and the normal turbulence intensity is decreased by the 
rotation. The experimental data show that the radial intensity is higher than 
the streamwise intensity. This relative magnitude is correctly predicted, 
although there exists a substantial difference in the magnitude of the individual 
stress component between the prediction and the experimental data. 

The predicted static pressure distributions at the mid-radius are 
compared with the experimental data in Figure 22. As the static pressure 
measurement itself contains some amount of error, it is difficult to comment 
as to whether the deviations between the predictions and the experimental data 
arise from the improper turbulence closure modelling or experimental error. 

The typical distribution of the rate of turbulent kinetic energy dissipa- 
tion is shown in Figure 23. These results show that the magnitude the decay 
rate of the energy dissipation is the same order of that of the turbulent 
kinetic energy. This indicates that the characteristic length scale of the 
turbulent flow field in the rotor wake (1 ^ k^''^/e) increases more rapidly 
than the rate at which the turbulent kinetic energy decreases. 


Comparison of Predictions with Two Different 
Turbulence Closure Models 


The rotor wake was also calculated with the so-called two-equation turbu- 
lence closure model (Launder and Spalding [15]) and the results were compared 
with those of the present turbulence closure model [modified Reynolds stress 
model, equations (39) through (44), (53), and (54)]. 





Figure 18. Profiles of Mean Velocity Components at 
Ravindranath and Lakshminarayana* [32]) 
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Figure 19(a). Reynolds Stress at r/rt*0.80 and s/C*»0.036 
(Data of Ravlndranath and Lakshmlnarayana , 
[32.3]) 






(Legend defined in Figure 18) 


o 



- 0.1 


(a) . Reynolds Stress a 
of Ravindranath ai 







70 


6 r 



Figure 21(a). Reynolds Stress at r/r^.^0.80 and s/C=0.036 (Data 
of Ravlndranath and Lakshmlnarayana, [32,3]) 
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Figure 22. Static Pressure Distribution at r/r^"0.80 (Data 
of Ravlndranath and Lakshmlnarayana, 132]) 
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Figure 23 . Energy 
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With the two-equation turbulence closure model, the local Reynolds 
stresses are evaluated with two local characteristic quantities of turbulence — 
the turbulent kinetic energy and the rate of turbulent kinetic energy dissipa- 
tion. In addition, two transport equations are provided to give local values 
of these two quantities. 

The two transport equations for the turbulent kinetic energy and the rate 
of turbulent kinetic energy dissipation in tensor form are as follows: 


where 
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and S is the source term in the transport equation of the rate of turbulent 
kinetic energy dissipation given in equation (55). 


In the present curvilinear coordinate. 
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dn ds 


The experimental data by Ravindranath and Lakshminarayana 132,3] were used for 
the comparison of results by two different turbulence closure models. The 
comparison of mean velocity is given in Figure 18. The agreement between the 
theory and the experimental data is good for both models; however, the modified 
Reynolds stress model predicts better results than the two-equation model. 

It should be mentioned here that the conventional k-£ model will not 
predict the shear stress or intensity components satisfactorily. The modified 
Reynolds stress model, however, can accurately provide this information. The 
comparison of Reynolds stresses are given in Figures 19, 20, and 21, where 
the present turbulence closure model shows substantial improvements in predic- 
tion of mean velocity as well as turbulent quantity. 

As explained earlier, the effect of rotation is to reduce the streamwise 
intensity and increase the radial component of turbulence intensity. The 
present model, which includes this effect, correctly predicts the trend (see 
Figures 19, 20, and 21). The radial component of intensities predicted by the 
k-e model is much lower than the measured data (Figure 20), while the author’s 
prediction comes closer to the data. 

The k-e model predicts higher values of shear stresses in the pressure 
side of the wake (Figures 19, 20, and 21). The mean strain rates are very 
high in this side, resulting in higher shear stresses when the k-e model is 
used for turbulence closure. 

These comparisons seem to indicate that the present Reynolds stress 
closure scheme models the rotation effect correctly through the rotation- 
originated redistribution term in the Reynolds stress transport equations. 

As explained in an earlier section, the rotor wake is highly three- 
dimensional and has multiple components of curvature. The transport equation 
of the rate of turbulent kinetic energy dissipation with the conventional source 
term, equation (55) is known to be responsible for the poor performance in 
predicting curved turbulent flows (Pope and Whitelaw [54], and Bradshaw [22]). 
This poor performance is reportedly caused due to inclusion of the mean shear 
rate in the production term. The energy dissipation equation with modified 
source term, equation (56) can predict curved turbulent flows better than that 
with the conventional source term. This point will be discussed further in 
the next chapter. 
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Comparisons of the predictions with two forms of the energy dissipation 
equation are made for the rotor wake data by Ravindranath and Lakshminarayana 
[32,3]. The comparison of mean velocity is given in Figure 24, Both turbulence 
closure models predict close values in mean velocity. The comparison of 
Re 3 rnolds stress is given in Figures 25, 26, and 27. Though there seems to be 
some improvement by the modified form of the energy dissipation equation, the 
difference is not substantial. 

Most of the curvature effect is sensed by the fluid particle near the wake 
center line where the velocity gradient is very large. The effect of curvature 
is directly related to the variation of angular momentum in the radial direction 
of curvature. For the rotor wake, the variation of angular momentum in the 
direction of the machine radius is small corpared to the variations in the 
direction of rg and Tg (Figure 2) , except for the hub wall and the annulus 
wall region. Furthermore the effect of curvature due to rg is opposite to that 
due to r^ near the wake center line (the sign of angular momentum change in 
equation (4) is opposite for rg and rg) . Therefore, the two curvature effects 
are cancelled and the total effect of curvature may not be substantial. The 
preceding argument holds only for the present case of rotor geometry and 
operating condition. Generally, the effect of curvature is significant. The 
decay of strearawise and radial velocity components predicted by the present 
modified Reynolds stress model with modified source term is compared with 
experimental data in Figure 28, where Uo, is the strearawise velocity component 
in the free stream, and U^, and are strearawise and radial velocity components 

at the center of the wake. Also, the decay of Reynolds stresses is compared 
in Figures 29 and 30. Where subscript m indicates the maximum value in the 
wake at that strearawise station. The agreement between predictions and experi- 
mental data is excellent, though some minor deviations are observed at the very 
near wake region. These comparisons indicate that the turbulence closure 
model employed in this thesis correctly models the rotation and curvature 
effects, through the modelled Reynolds stress transport equation. 
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Figure 27. Reynolds Stress at r/r|-«0.80 and s/C*0.036 (Data 
of Ravlndranath and Lakshinlnarayana, [32,3]) 






82 



(J^/U)j^ x1Q2 




Figure 29. Decay of Maximum Turbulence Intensity (Data of 
Ravindranath and Lakshminarayana, [32,3]) 
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Figure 30. Decay of Maximum Shear Stress (Data of Ravindranath and 
Lakshmlnarayana , [32 , 3 ] ) 
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NUMERICAL STUDY OF ASYMMETRIC TURBULENT WAKES 
OF SINGLE AIRFOILS AND CASCADES OF AIRFOILS 


Previous Analytical and Experim^tal Investigation 
of the Isolated Airfoil aiid Cascade Wake 


As free shear flows, the wakes of single airfoils and of cascades of air- 
foils have their own physical natures. The mean velocity profiles of the single 
airfoil wake and the cascade wake is two-dimensional and generally asymmetric. 

The asymmetric nature in the wake is due to the past history of the flow. For 
most engineering applications, the Reynolds number, based on the chord length 
of the airfoil, is large enough for the boundary layer near the trailing edge 
to be considered fully turbulent. Therefore, the subsequent wake is also 
turbulent. The cascade wake has a periodic distribution of flow quantities 
with a period equal to the spacing of the blades. Both the single-airfoil wake 
and the cascade wake have an elliptic nature If the flow field very near to the 
trailing edge of the airfoil is considered. The static pressure change in the 
single airfoil wake has been reported to be mild both in normal and streaniwise 
direction (Chevray and Kovasznay [55], and Fermin and Cook [56]). For the 
cascade wake, the static pressure is considered to have an adverse gradient 
in the streaniwise direction because the edge velocity in the cascade wake 
decreases downstream. 

The understanding of the turbulent wakes of the single airfoil and the 
cascade of airfoils represents not only an interesting problem in fluid mechanics, 
but also an important aerodynamic problem related to the prediction of lift, 
drag, and aerodynamic noise level of the airfoil and the cascade. 

Very few experimental and theoretical studies have been made for the 
turbulent wakes of the single airfoil and the cascade of airfoils. Near and 
far wakes of a symmetrical airfoil were first investigated experimentally by 
Silvers tein et al. [57], who provided empirical relationships for the wake 
decay. Preston et al. [58] carried out a systematic investigation of the 
characteristics of the wake behind an isolated airfoil and observed that a 
similarity in mean velocity profile exists close behind the airfoil. 

Mendelsohn [59] measured the wake of an isolated airfoil at various incidence 
angles. All of these measurements were concentrated on the mean velocity, and 
turbulence flow quantities were not measured. Chevray and Kovasznay [55] 
measured the mean velocity and turbulent quantities in the symmetric wake of 
a flat plate including the trailing edge region. Though their experimental 
data have been frequently referred to as the most accurate and comprehensive 
data for the two-dimensional turbulent wake by many investigators (Launder, 

Reece and Rodi [30], Bradshaw [60], etc.), the experimental data do not 
include the measurement of spanwise intensity. No turbulence measurements for 
the asymmetric wake of a single airfoil including the near wake region have 
been reported before this report. The turbulent structure as well as the mean 
velocity profile of an asymmetric turbulent wake is known to be different from 
that of a symmetric wake. Measurements of turbulence quantities as well as 
the mean velocity in an asymmetric wake of a single airfoil will be reported in 
the following sections. 
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Theoretically, most work related to the two-dimensional turbulent wakes 
has been confined to the far wake region. This is because the flow evolution 
far downstream is very slow and a subsequent similarity solution can be assumed 
(Tennekes and Lumley [4J). Though the similarity solution has good agreement 
with experimental measurements, it cannot be easily expanded to the upstream 
flow. The measurements in the turbulent wake of a cylinder by Townsend [61] 
shows the similarity in mean velocity about 80 cylinder diameters behind and 
the similarity in Re 5 molds stress beyond about 200 cylinder diameters. Further- 
more, the more complicated the statistical qiiantity, the longer it takes to 

reach the region which is characterized by the similarity rule. For the near- 

wake region, Goldstein [5] obtained a semi-analytical solution for the 
symmetrical laminar wake of a flat plate. Also, the symmetric laminar near 
wake of a flat plate including the trailing edge region has been successfully 
predicted with various numerical schemes (Dennis and Dunwoody [62], Dennis and 
Chang [63], and Stewartson [6]). Compared to the laminar wake, the turbulent 

wake has a much more complicated structure at the trailing edge region. The 

boundary layer at the trailing edge of an airfoil has a viscous sublayer, 
buffer layer, inertial sublayer, and an outer layer. The transfer from the 
boundary layer structure at the trailing edge to the wake structure occurs in 
a short distance and is not well-understood. 

Numerically, the symmetric turbulent wake of a flat plate has been 
predicted by Launder, Reece, and Rodl [20] and Pope and Whitelaw [54] with 
slightly different numerical schemes. As the asjTiimetric turbulent wake has 
a different turbulence structure from that of the symmetric wakes, the turbu- 
lence closure modelling for the symmetric turbulent wake is not considered to 
be directly applied for the predictions of asymmetric turbulent wakes. Proper 
modification of existing turbulence closure models for the prediction of 
asymmetric turbulent wakes will be discussed and the numerical results will 
be compared with the experimental study in this thesis. 


Numerical Analysis of Turbulent Wakes of 
Isolated Airfoils and Cascades of Airfoils 


The Curvilinear Coordinate System 

As was evident from the previous flow visualization around a single air- 
foil, there exists streamline curvature in the near wake when the flow field is 
not symmetric. For calculation of the curved turbulent wakes of a single air- 
foil and a cascade, the curvilinear coordinate system shown in Figure 31 was 
introduced. In this figure, r is the distance from the origin of curvature, 
s is the strearawise distance on the specified path, and b is normal to s and r. 
With this coordinate system, the actual streamline of the wake can be closely 
represented. This coordinate system forms part of the coordinate system for 
the calculation of the wakes of the single airfoil and the cascade blade. The 
fundamental metric tensor, Christoff el symbols can be derived in the same way 
as for the curvilinear coordinate system for the rotor wake. Detailed deriva- 
tions are not repeated here. The fundamental metric tensor for this coordinate 
system Is: 
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Chris toff el symbols of the second kind are as follows: 
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All other Chris toff el symbols of the second kind are zero. 

The Gove rning Equation in the Curvilinear Coordinate System 

The equations governing the steady incompressible flow are introduced in 
generalized tensor form. 
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The continuity and momentum equations are, respectively, 
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In the present coordinate system, the continuity equation is: 
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Where U, W are mean contravariant velocity components in the present curvi- 
linear coordinate system and u, w are corresponding fluctuating velocity 
components. Two components of the equation of mean momentum conservation in 
the present curvilinear coordinate system, at high Reynolds number are: 
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Turbulence Closure Modelling for the Curved Turbulent Flows 

The turbulent wakes of a single airfoil and a cascade which develop under 
the influence of streamline curvature were calculated using three different 
turbulence closure models. The first model was comprised of transport equations 
for the turbulent kinetic energy and the rate of energy dissipation. The second 
model used equations for the rate of turbulent kinetic energy dissipation and 
Reynolds stresses, but the effects of the convection and diffusion in the 
Reynolds stress transport equation were handled collectively. The third model 
utilized equations of turbulent kinetic energy dissipation and Reynolds 
stresses in nearly exact form. All three of the models were modified for the 
effects of streamline curvature. 
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The transport equations for the turbulence closure which, make up the 
three turbulence closure models are. as follows: 


First model: two-equation model. 
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and S is the source term in the transport equation of the rate of turbulent 
kinetic energy dissipation. Two forms of S were utilized and the results by 
each form were analyzed. The two forms of S are: 
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Second and third models: Reynolds stress models. 

The second and third models utilize the simplified form of the transport 
equation of Rejmolds stress. The pressure-strain correlation term in the 
transport equation of Reynolds stress was modelled following the proposals by 
Rotta [50] and Naot, Savit, and Wolfstein [51]. Furthermore, the rate of 
turbulence energy production, which appears in the modelling of mean-strain 
effects in the pressure-strain correlation, was replaced with the rate of 
turbulent kinetic energy dissipation. The second and third models differ only 
in the handling of convection and diffusion terms. The convection and 
diffusion terms are handled collectively in the second model; the effects of 
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the two terms are assumed to be proportional to the production term and the 
magnitude of the effects are corrected at each, iteration. 

The third model includes the individual effects of convection and diffusion 
and is similar to the equation used by Launder, Reece, and Rodi [20] for 
straight non-rotating flows. The rotation-originated redistribution term was 
added in the turbulence closure equations. 


Second model: simplified Reynolds stress model. The equation used for 

this case is given as follows; 
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where relates the collective effects of the convection and diffusion terms 
to the production terms. 


Third model; full Re 3 molds stress model. The complete Reynolds stress 
equation (with suitable modelling for pressure-velocity correlation, etc.) 
is given as follows ; 
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The values of constants in equations (76) through (81) were not optimized for 
the present curved wakes; the values used by Pope and T-Thitelaw [54] were 
adopted. The value of in equation (79) was 1.8 for the present calculation. 

The effect of streamline curvature on a fluid flow has been qualitatively 
explained with the motion of a disturbed element of fluid. The argument based 
on Incompressible, inviscid fluid is as follows. The effect of streamline 
curvature is to stabilize the flow field and consequently diminish the turbulent 
transportation when the angular momentum increases with the radius of streamline 
curvature and vice versa for the opposite case. Quantitatively, the following 
expression can be derived for the criterion of stability: 
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where O) is the corresponding frequency, U is the circumferential velocity 
component, and r is the radius of curvature. For most thin shear flow, the 
first term is dominant. To account for the above effect of curvature, Bradshaw 
[22] proposed modification of the conventional effective viscosity model and 
Launder, Pridden, and Sharma [24] modified the two-equation model. Irwin and 
Smith [64] reported successful calculation by manipulating curvature terms which 
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appear In the transport equation of Reynolds stress In the special coordinate 
system. The turbulence closure models with the conventional dissipation equa- 
tion (with the source term (81)) have been widely applied and criticized for 
poor prediction of curved flows Ce«R. t Launder, Vridden, and Sharma, [24], and 
Bradshaw [22]). 

The existing prediction techniques for curved turbulent flow utilizing 
the turbulence models with the conventional dissipation eqiiation (with the 
source term (78)) are considered to be inadequate. This is caused by the 
inappropriate form of the dissipation equation which is based on unproven 
assumptions and tested with only straight flows. 

To represent the effect of streamline curvature properly, the rate of 
energy dissipation should be increased where the angular momentum increases with 
the radius of curvature and decreased where the angular momentum decreases with 
the radius of curvature. For the coordinate system shown in Figure 31, the 
production term in equation (76) is, neglecting high-order stress terms. 
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and if higher order curvature terms including s/r are neglected, 
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and for the first model, 
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The above expression can be further simplified for the thin shear layer with 
small pressure gradient in the normal direction. 
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The second term of the above expression has the same form as the principal 
term of expression (82) and represents the curvature effect that was discussed 
qualitatively. As will be shown in the next section, turbulence closure models 
with the conventional dissipation equation [with the source term (78)] do not 
predict curved turbulent wakes accurately even though the transport equations 
are in the curvilinear coordinate system. Referring to equation (76), it is 



92 


clear that the turbulent kinetic energy increases with a decrease in angular 
momentum CcU). However, the dissipation rate also Increases when the production 
of kinetic energy dissipation rate, and the overall effect of curvature is not 
adequately Included in the turbulence modelling. 

Chambers and Wilcox [25], proposed modification of the first model for the 
curvature effect. They suggested adding extra singular terms to the transport 
equations of turbulent klnetit energy and the rate of energy dissipation. 
However, these terms increase both the kinetic energy and the rate of energy 
dissipation where the angular momentum decreases with the radius of curvature, 
and subsequently do not fully account for the curvature effect for the present 
set of data. 

Two different dissipation equations were tried for better prediction of 
curved turbulent wakes. As Lumley and Khajeh-Nouri [21] argued, the production 
term can be used as the source term in the dissipation equation with the 
assiimption that 


(uV - I . (85) 

This assumption is not valid when the flow field is affected by streamline 
curvature and rotation; therefore, the exact source term, 


(uV - 4 k6. .Ku^u"^ - -| k6..) 

by Lumley and Khajeh-Nouri [21] was utilized. To make the dimensions correct, 
the following source term was used for the present calculation; 


S = C. - I k6^.) (u^u^ - I k6 ) (86) 

which is given as the second form of the source term in equation (77) . 

Also, a modification of the conventional dissipation equation was tried 
in order to represent the curvature effect more accurately. The coefficient 
of the source term, which includes production of energy, was modified ‘as 
follows: 
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where R . = — — and can be interpreted as the Richardson number in curved 
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relates the streamline curvature to the local 


turbulent flow. The constant 
dissipation rate, and the optimized value with the available data was 0.03. 
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Nine turbulence closure models used for calculation of the wakes of 
Isolated airfoils and cascades of airfoils are tabulated and shown in Table 2. 


Numerical Technique 

The niunerlcal scheme for the turbulent wakes of single airfoils and 
cascades of airfoils is very similar to that used for calculation of the rotor 
wake. The wakes of single airfoils and cascades of airfoils are two-dimensional, 
and the numerical scheme developed for the rotor wake can be simplified for the 
calculation of two-dimensional wakes. The continuity and mean momentum equa- 
tions (73), (74), and (75), along with the closure equations (76) through (81) 
were solved in elliptic form. At the near-wake region, the streamwlse velocity 
gradient is of the same order of magnitude as the normal velocity gradient. 

None of the shear stress conqjonents are negligible. Therefore, the conventional 
boundary-layer assumption cannot be invoked in this region. The parabolic 
marching technique does not produce accurate predictions. The initial values 
of unknowns were obtained by parabolic marching technique, then the values 
were corrected iteratively. 

The t 3 Tpical domain of calculation of the wake of an Isolated airfoil is 
shown in Figure 32. The boundary conditions on AB and CD were based on the 
experimental data. The turbulent kinetic energy on boundaries AD and BC was 
calculated on the basis of the free-stream turbulence intensity (0.2%). The 
boundary values of dependent variables except the rate of the turbulent kinetic 
energy dissipation were directly obtained from the experimental data. As 
mentioned earlier, the rate of turbulent kinetic energy dissipation was assumed 
to be equal to the production of the energy. For the first iteration, the 
origin of the coordinate system in Figure 31 was located on the line which is 
normal to the camber line and passes the trailing edge of the blade* This 
distance between the origin of the coordinate system and the trailing edge was 
assumed to be one chord length for the first calculation. The iteration was 
performed downstreamwise from the trailing edge. After each sweep of the 
Iteration, the origin of the coordinate system was corrected and the effect 
of the streamline curvature was more accurately estimated for the next itera- 
tion. The turbulent kinetic energy and rate of energy dissipation were handled 
as diffusion terms in the first and second model. In the second model, the 
Reynolds stress components were first computed using eqtiation (75) and the 
values of U^, k, e at the previous iteration. The governing equations were 
solved with these Reynolds stress components. Typical positions of grid nodes 
are given in Figure 33. In the third model, the shear-stress nodes were 
located between the adjacent nodes in the normal direction to provide stability 
in calculation. About 20 iterations gave converged values for all closure 
models. Figure 34 shows calculated values at different iterations for the wake 
of a single airfoil with the second turbulence closure model. 

About 20 seconds were necessary for each Iteration with the first model on 
IBM 370/3066 at The Pennsylvania State University. The second and third models 
required about 20 and 50 percent more computing time, respectively. 




Turbulence Models 


Table 2. Turbulence Closure Models with Different Forms of the Source Term in the Equation 
for the Rate of Turbulent Kinetic Energy Dissipation 
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A well-documented problem, the turbulent wake of a flat plate was solved 
to test the present numerical scheme. This symmetric turbulent wake, measured 
by Chevray and Kovasznay 155] was previously calculated by Launder et al. [20] 
and Pope and Whitelaw [54]. The predicted variation of the wake centerline 
velocity is shown in Figure 35, along with the experimental data. The three 
turbulence closure models predicted almost identical mean velocity profiles. 
Launder et al. [20] and Pope and Whitelaw [54] also reported almost identical 
results for this standard problem. Therefore, no additional numerical results 
for this problem are compared here. 


Comparison of the Numerical Predictions of the Isolated 
Airfoil Wakes with Corresponding Experimental Data 


The turbulent wakes of the isolated airfoil were numerically calculated 
and the results were compared with the experimental data. The three different 
turbulence closure models, equations (77) and (78), (80), and (81) were used 
for the calculation, along with three different forms of the source terms in 
the dissipation equations, equations (78), (86), and (87). Hence, the numerical 
predictions with nine different sets of turbulence schemes were compared with 
the experimental data. The turbulence closure models are shown in Table 2. 

The experimental data by Hah and Lakshmlnarayana [52] were used for the 
comparison. The measuring stations for the single airfoil wake are shown in 
Figure 35. 

The calculated profiles of mean velocity and shear stress without curva- 
ture modification are given in Figure 36. The experimental data are also 
provided for comparison. The results with curvature modification are presented in 
Figures 37 and 38. The figures show that the three turbulence models do not predict 
the mean velocity and shear stress profiles accurately when the dissipation 
equation with the conventional source term (78) is used. All three models 
predict higher shear stress at the suction side than was foiond from the 
experimental data. However, this discrepancy disappears by modifying the 
conventional source term, equation (87) or by adopting the exact source term, 
equation (86). All three models predicted almost identical results for the mean 
velocity and the shear stress when the dissipation equation was modified. 

The calculated values of turbulence intensity with different closure 
models are shown in Figures 39 through 47. The first model predicts almost 
isotropic values for the three intensity components, and the profile is corrected 
when the modified or exact source term is used in the dissipation equation. 

The second and third turbulence closure models with the conventional source term 
in the dissipation equation predict a higher intensity level on the pressure 
side and lower intensity level on the suction side than was found in the 
experimental data. However, the discretancy between prediction and experiment 
is considerably decreased by modifying the coefficient of the conventional 
source term or by adopting the exact source term. As described earlier, the 
turbulence level in the wake is not linearly dependent on the mean shear rate. 

The turbulence level on the pressure side is very low compared to that on the 
suction side of the wake, even though the mean shear rate on the pressure side 
is higher than that on the suction side. When the production of turbulent 
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Figure 39. Relative Streamwise Intensity at s/C==0.032 
and i=6® 
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Figure 40. Relative Streamwise 
and i-6P 
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Intensity at s/C=0.032 
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Figure 41. Relative Streamwise Intensity at s/C*0.032 
and i=*6® 
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Figure 43. Relative Normal Intensity at s/C*i 
and 1*6® 
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Figure 45. Relative Spanwise Intensity at s/C-0.032 
and i“6® 
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kinetic energy Is used as the source term In the dissipation equation, the 
turbulence quantities become strongly linked to the mean-shear rate and the 
prediction becomes Inaccurate. When the modified or exact source term are 
adopted in the dissipation equation, this strong link between the mean shear 
rate and the turbulence quantities no longer exists, and more accurate predic- 
tions are obtained. The predicted overall decay rate of mean velocity is shown 
In Figure 47b. 


Cascade Wake Data and Comparison with Predictions 


The wakes of a cascade measured by Raj and Lakshminarayana [53] were 
calculated for further comparison. The blade profiles were similar to the 
NACA-65 ( 8 A 2 section, and other important parameters are: 

solidity C/S = 1.505, inlet angle = ^5®, incidence 1 - 2®, 0®, - 6 ®, 

and Reynolds number based on chord length - 3 x 10^. 

The upstream velocity and free stream turbulence intensity were 24 ms"^ and 
0.16 percent. Because the experimental data did not provide spanwise inten- 
sity, the energy distribution at the inlet boundary was guessed with available 
single airfoil wake data. Also, the rate of turbulent kinetic energy dissipa- 
tion was assumed to be equal to the production of energy. Because of the 
characteristics of the flow around a cascade, the periodic boundary condition 
was applied between corresponding points in normal direction. The flow field 
and coordinate system Is shown in Figure 48. 

Since the experimental data of i = -6 provides the most detailed informa- 
tion on the boundary conditions, the comparison is made for the experimental 
data of i = - 6 ®. Figures 49, 50, and 51 show predicted mean velocity and 
Reynolds stress distribution with three turbulence models. The modified 
turbulence closure models show substantial Improvement in prediction compared 
to original closure models, and all the three models predict almost identical 
results if they are properly modified. 

The predicted intensity distributions are shown in Figures 52 through 57. 
The first model predicted only isotropic intensity values, though the trend 
was improved by modifying the conventional source term for curvature or adopting 
the exact source term in the dissipation equation. The second and third models 
provided better predictions for individual components, and Improvements were 
observed when the dissipation equation was modified. 

The discrepancy between theory and experiment in the distribution of 
streamwise intensity near the center of the wake is not explained at this 
stage. It may be due to scatter in measurements or to improper turbulence 
modelling. The trajectory of the wake centerline is shown compared with 
predicted values in Figure 58. Even though the curvature seems moderate, its 
effect on turbulence structure, as well as mean velocity, is appreciable. 
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Figure 49. Streamwlse Velocity and Shear Stress Profile 
at s/C*0.08 and l“-6** (Data of Raj and 
Lakshmlnarayana, [53]) 










Figure 50 , Streamwlse Velocity 
at s/C=0.08 and i=-( 
Lakshmlnar ayana , [5 
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Figure 52. Relative Streamwlae Intensity at s/C"0.08 and 

(Data of Raj and Lakshmlnarayana , [33]) 
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Figure 53. Relative Streamwlse Intensity at s/C*0.08 

and i«-6® (Data of Raj and Lakshmlnarayana, [53]) 
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0 Expt . (Ref. 53) 

Prediction with Model A3 

Prediction with Model B3 
Prediction with Model C3 



Figure 54. 


Relative Streamwise Intensity at s/C=0.08 

and i=“6® (Data of Raj and Lakshminarayana, [53]) 
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Figure 57. Relative Normal Intensity at s/C«0.08 

and 1*-^ (Data of Raj and Lakshmlnarayana, [53]) 
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FORTRAN PROGRAM FOR THE PREDICTION 
OR rotor wake DEVELOPMENT 


The Fortran computer program that calculates the development of the turbu- 
lent rotor wakes In a compressor Is designed to be run on the IBM 370 computer 
system, but should prove easily adaptable to other installations. It is 
approximately 1450 cards in length and consists of a main program and one 
short subroutine. The subroutine ELIM is to solve tri~diagonal matrices and 
any efficient routine can be used instead of this routine. Various routines 
to generate curve-fitted data of experimental results are not included in 
this appendix because any kind of utility subroutines for this kind of mathemat- 
ical operations can be used. Various subroutines used for plotting the results 
are not included for the same reason. Some portion of computing time and 
computer storage may be saved by several modifications in the computer code, 
but this has not yet been tried. 

The development of turbulent wakes of a single airfoil and a cascade of 
airfoils can be also predicted with the present computer program. The wake of 
a cascade of airfoils can be considered as a simplified case of rotor wake, when 
the radial distance is very large and the rotational speed is zero. The 
single airfoil wake is a further simplification of the rotor wake, with large 
blade spacing or very small chord length compared to the blade spacing. There- 
fore the present computer program can be used for the prediction of a cascade 
and an airfoil wake. 


Main Program 


The main program first reads input data cards, which are detailed later 
in this section. As the present numerical scheme is for elliptic calculation, 
boundary conditions are required for all the boundaries of the calculating 
domain. The input data includes all these boundary conditions. The geometry 
of the compressor, the angle 8 at various radial stations, and other constants 
are given as input data. The present program is written to receive boundary 
conditions as a combination of polynomials and Gaussian functions. However, 
this is quite arbitrary and any kind of curve-fitted distribution of unknowns 
on the boundary surface can be used. Seven equally spaced radial stations are 
used for the finite difference calculation. The normal stations vary from 75 
at first radial station to 87 at the last radial station. At each radial 
station, two different spacings between normal station are used. In the 20 
nodal points from the central nodal points in both normal directions, the 
spacing between nodal points are 1/4 of the spacings of outer radial points. 
Also the present calculation domain contains one complete rotor passage with 
the blade trailing edge at the center. Therefore, the normal spacing varies 
in the radial direction to cover the spacing between blades with the given 
number of radial points. The number of radial stations and normal stations 
in finite difference calculation is also arbitrary even though it requires some 
change of constant in the program and can be done without any major effort. 

Any type of automatic or numerical mesh generation in the surface normal to the 
strearawise direction can be incorporated in the present computer program 
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without any difficulty. With the given Initial and boundary conditions, the 
six components of Re3molds stress are calculated and with this Reynolds stress 
components, the momentum and turbulence closure equations are solved. The 
detailed operations in the main program are shown in the flow chart (Figure 6) . 
Because the available computer memory storage was quite small (around 200 k 
with IBM 370/3033 at The Pennsylvania State University), the computed results 
were stored in the RJE batch files and the unknown quantities stored in the 
batch files were updated during numerical calculation. 


Input Data Cards 


Input data cards are as follows. Various formats are also given for 
each data card. 


Card 1 Card 1 contains one integer and two floating point numbers 
(15, 2F10.4) 


Columns 

Name 

Description 

1 

to 

5 

lAL 

Number of iteration 

6 

to 

15 

DENS 

Density of fluid (kg/m^) 

16 

to 

25 

SOR 

Successive over-relaxation factor. 


A fixed value of 1.55 was used 
for the present study. 


Card 2 Card 2 contains three floating point numbers (3F10.3) 

Columns Name Description 

1 to 10 DR Radial spacing in finite difference 

grid (m) 

11 to 20 DS Streamwise spacing in finite difference 

grid (m) 

21 to 30 DN Normal spacing in finite difference 

grid (m) 

Card 3 Card 3 contains two floating point numbers (3F10.3) 


Columns 

Name 

Description 


1 

to 

10 

EF 

Turbulence kinetic energy at free 
stream (m^/sec^) 

11 

to 

20 

DR 

The rate of turbulence 
dissipation at free 

kinetic energy 
stream (m^/sec 


Card 4 Card 4 contains one floating point number in columns 1 to 10, 
the angular speed of rotor. 


Card 5 


Card 5 contains one floating point number in columns 1 to 10, 
AL, the radial variation of the angle between sterarawise direc- 
tion and machine axis (radian/m) . 
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Card 6 


Card 6 contains four floating point numbers for the descrlp' 
tlon of the compressor geometry (4F10.3) 


Columns Name 


Description 


1 to 10 VNON(l) 
11 to 20 B(l) 

21 to 30 SA(1) 

31 to 40 R(l) 


Free-stream velocity at hub (m/sec) 
Outlet angle at hub (Radian) 

Blade spacing at hub (m) 

Hub radius (m) 


Card 7 


Card 7 contains three floating point numbers to describe 
radial variation of unknown quantities. 


Columns Name Description 

1 to 10 BCl Radial variation of outlet angle (radian) 

11 to 20 SACl Radial variation of spacing (m) 

21 to 30 VNONCl Radial variation of free-stream 

velocity (m/sec) 


Cards 8 
to 18 


Cards 8 to 18 are used to give coefficients of functions 
for the distribution of the streamwise velocity on the 
boundary surfaces. The following function was used for the 
description of streamwise velocity component 


U(n) = A^U - A«U 

' ' 1 oo 2 ° 


“ -0.693n/Lg 
e ® - 


1.0 




+ A,(^)^ 




A 


1 


where Uqo = free stream velocity 

Lg = length scale, the distance from the wake center to 
the point where the velocity defect is half of 
that at the wake center. Different length 
scales are used in the pressure and suction 
sides of the wake. Lg = A 3 in suction side and 
Lg = A^ in pressure side, 
n = normal distance, n = 0 , at the wake center 
,A 2 *’*’ Ay = coefficient of functions 


Each card contains seven floating point numbers to describe 
the seven coefficients used. 


Columns 

Name 

Description 

1 

to 

7 

xcd.i.i) 

Al 

8 

to 

14 

XC(1,2,I) 

Ai 

15 

to 

21 

XC(1,3,I) 

A 3 

22 

to 

28 

XC(1,4,I) 

A 4 

29 

to 

35 

XC(1,5,I) 

A5 

36 

to 

42 

XC(1,6,I) 

H 

43 

to 

49 

XC(1,7,I) 
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Cards 19 Cards 19 to 29 are used to give coefficients of functions 
to 29 for the distribution of the radial velocity on the boundary 

surfaces. The following functions were used for the descrip- 
tion of the radial velocity component. 

W(n)/Uoo - WD4 if n/S < X34 

W(n)/U^ » (WD4 - WD3)x(n/S)/(x3^ - X33) + WD4 

- (WD3 - WD4)x(x33 ^^^*34 “ ^33^ 

if = n/S < +X33 

W(n)/U^ = +(WD3 - WC)x(n/S)/x33 + WC if +X33 ^ n/S < 0 

W(n)/U^ = (WD^ - WC)x(n/S)/x3^ + WC if 0 ^ n/S < 

W(n)/U^ = (WD^ - WD^)x(n/S)/(x32 “ + WD2 

-(WD 2 - WD^)x(x32)/(x32 " ^ 3 ;!^) 
if < n/S < X32 

W(n)/U^ “ WD^ if n/S > x^^ 

where W(n) = radial velocity 

S * spacing of the blade 
n - normal distance from the wake center 
WC,WDi,WD 2 , 

WD3,WD4,X3i, 

*32»^33»^34 * constants 

Each card contains nine floating point numbers to describe 
nine coefficients used. 


Columns 

Name 

Description 

1 to 6 

XC(3,1,I) 

WC 

7 to 12 

XC(3.2,I) 

X3I 

13 to 18 

XC(3,3,I) 

TOx 

19 to 24 

XC(3,4,I) 

X32 

25 to 30 

XC(3,5,I) 

WD2 

31 to 36 

XC(3,6,I) 

X33 

37 to 42 

XC(3,7,I) 

WDo 

43 to 48 

XC(3,8,I) 

^34 

49 to 54 

XC(3,9,I) 

WD4 
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Cards 30 Cards 30 to 40 are used to give coefficients of functions for 
to 40 the description of the normal velocity on the boundary 

surfaces. The following functions were used for the descrip- 
tion of the normal velocity component. 

V(n)/U^ " if < +*24 

V(n)/U„ - +(V^ - V3)x(n/S)/(X24 - *23) ^^4 

-(V^ - V3)x(x24)/(X24 " *23^ "*24 ^ ^ '^*23 

V(n)/U^ * +V^ X (n/S)/x23 If +x^3 « n/S < 0 
V(n)/U^ * X (n/S)/x2i If 0 < n/S < 

V(n)/U^ = (V^ - V^)x(n/S)/(x22 “ ^ 21 ^ + ^2 

-(V^ - V^)x(x22)/(x22 ^ *21^ 

If *22 

V(n)/U^ - V2 if X22 ~ n/S 

where Ucx> = free stream velocity 
S = spacing of the blade 
n = normal distance from the wake center 
V = normal velocity 

Vi,V2,V3.V4, 

^21»^22»*23» 

X24 = constants 

Each card contains eight floating point numbers to describe 
eight coefficients used. 


Columns 

Name 

Description 

1 

to 

6 

XC(2,1,I) 

*21 

7 

to 

12 

XC(2,2,I) 

Vl 

13 

to 

18 

XC(2,3,I) 

*22 

19 

to 

24 

XC(2,4,I) 


25 

to 

30 

XC(2,5,I) 

*23 

31 

to 

36 

XC(2,6,I) 

V3 

37 

to 

42 

XC(2,7,I) 

*24 

43 

to 

68 

XC(2,8,I) 

V4 
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Cards 41 Cards 41 to 51 are used to give coefficients of functions 
to 51 for the description of the turbulence kinetic energy on the 
boundary surfaces. The following functions were used for 
the description of the turbulence kinetic energy. 

EN(n)/U^ = E4 if n/S < - 

EN(n)/U^ = (E4 - E3)x(n/S) / (-x^^ + x^^) + E4 
-(E4 - E3)x(-x^^)/(-x^^ + x^^) 

if ^ n/S < -x ^3 

EN(n)/U^ = (E3 - E5)x(n/S)/ (-x ^3 + + E3 

-(E3 - E5)(-x^2)/(-x^ 2 + x^^) 

If -x ^3 i n/S < -x ^5 

EN(n)/U^ = (El - E5)x(n/S)/(x + x^j) + El 

~(E1 - E5)x(x^^)/(x^j^ + x^^) 

if -X, c = ^/S < X, - 

45 41 

EN(n)/U^ = (E 2 - El)x(n/S) / (x ^2 “ ^ 4 ^) 

-(E2 - E1)x(x^2^/(x^ 2 “ ^41^ 
if x^^ < n/S < x ^2 

EN(n)/uf = E2 if x^^ = 

where Uqo = free stream velocity 

EN = turbulence kinetic energy 
n = normal distance from the wake center 
E1,E2,E3,E4,E5, 

^41»^42»^43* 

X 44 >X 45 *= constants 


Each card contains ten floating point numbers to describe 
ten coefficients used. 
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Columns 

Name 

Description 

1 

to 

6 

XC(4.1,I) 

X41 

7 

to 

12 

XC(4,2,I) 

El 

13 

to 

18 

XC(4,3,I) 

*42 

19 

to 

24 

XC<4,4,I) 

E2 

25 

to 

30 

XC(4,5,I) 

X45 

31 

to 

36 

XC(4,6,I) 

E5 

37 

to 

42 

XC(4,7,I) 

X43 

43 

to 

48 

XC(4,8,I) 

E3 

49 

to 

54 

XC(4,9,I) 

X44 

55 

to 

60 

XC(4,10,1) 

E4 


These cards are not utilized in the Fortran program listing attached. 
The program can be easily modified to include the generalization of the 
input data mentioned above. 

Because many experimental data indicate that the distribution of 
turbulence kinetic energy can be approximately represented with the mean 
shear rate when proper coefficients are used. This simple representation 
of turbulence kinetic energy is used instead of the above multi-functions 
for the same input data. 
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Typical Input Data Set 

R.01^ 0.0014 0,00£Elo 

2.0 4.0 
111.63 

1 . 286 

29.6 0,502 0.1018 0.341 

0.0225 0.0051 1.7 

0 . 209 0 . 745 0 . 003 0 . 00 1 ~0 . 00260 . 000 1 7-0 . 080 1 

0 . 0 
0.0 
0 . 0 
0 . 0 

0.392 0.481 0.003 0.004 -0.017 0.064 0,6003 

0,426 0.349 0.017 0.624 0,0186 -0.067 -0.0349 

Ef . 243 0 . 723 E' . 003 6 , 603 -Ei . i3 1 -0 . C'Ef0Ef0 . 00004 

0.632 0.312 0.003 0.005 -0. 6522-0. 00136. 60036 

0 . 593 0 . 352 0 . 622 E' . 628 6 .016 2 -6 . 62 1 ~ 0 . 6029 

0.573 0.411 0.026 0.04 -6.094 0.666 0.069 

-6 . 289© . E'62 6 . 387 6 .5 0 . 043 -6 . 0620 . 1 26 - E? .27 6 . 066 

0 . 6 
0.0 
0 . 0 
0 . 0 

-6.1270.094 6,333 6.5 0.049 -0.12 0.235 -0.3 6.065 

6,016 0.117 0.16 9 0 , 882 6 . 055 -0 . 1 946 . 234 -0.7146. 117 
-6.3826.074 6.56 1.0 0.164 -0.1 0.256 -0.5 0.155 

6.125 6.1 0.281 1.25 0,127 -6.16 0.3 -0.6 0,163 

-6.1870.2 0.096 1.0 0.023 -6,2 -6. H 1-^0, 5 -0,037 

-6.18 0.2 0.631 1,0 0.02 -0.5 0.062 -1.0 0.222 

6.042 -0.1270,1 -0,014-0.03 6.101 -0.33 6.015 

0 . 0 
0.0 
0.0 
e.e 

0.105 -6.13 0,5 -0.03 -6.1320.673 -0.55 0.02 

, 29 —0. 15 0. (■' — 0 , 022— 0 . 254E* . 0y —0. I'' 0r 02 

0.652 -6.J4 0.12 -0.055-6.0610.137 -0.4 0,1 

6.1 -0.1646.3 -6,68 -6.1 0.2 --0.5 0,22 

6,26 -0.1 0.72 -0.01 -0,03 0.078 -1.0 

0.21 -0.0870,74 -6.004-6.5250.076 -1.0 

0.64 0.12 0.13 -6.0? 0,011 0.03 0.04 0.06 -E.PJ -0.025 

0.0 
0.0 
6.0 
0.0 


0. EtAp: 

0, 

142 

6 . 

11 

-0.05 

Er . 05 

0.621 0.03? 

6 . 0?6 

-6.05 -0,001 

6 . 67 

0. 

15 

_ 

12 

~6,0E 

0.05 

0,02 0.08 

0.025 

-6.06 -6.061 

6.05 

0. 

12 

0, 

128 

-6.02 

0. 03 

0.002 0,06? 

0.0?4 

-6.02 -0.004 

6.05 


10 

0. 

13 

-0.64 

6.034 

0.002 0.05 

0.04 

-6,05 -6.004 

0.11 

El . 

15 

0, 

11 

-6.3 

0,03 

0.00150. 11 

0 . 0 1 87 

-6.? -6.001 

0. 1 1 

E:". 

12 

0. 

11 

-6. 3 

0.03 

6.00150. 11 

0.0187 

-6.3 -6,601 
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CONCLUSIONS 


The major conclusions derived from this Investigation are as follows: 

1. The development of the turbulent wake of rotor blades In a 
compressor was predicted numerically. The numerical results 
show excellent agreement with various wake data from different 
blade geometries. All the components of the mean velocity profile, 
as well as the turbiilence Intensity and stress profiles, are 
predicted accurately. 

2. The non-orthogonal curvilinear coordinate system utilized describes 
the actual streamline closely, and provides simplicity In describing 
the governing equations and the boundary conditions. 

3. The Reynolds stress transport equation was modified for the 
effects of rotation and the streamline curvature. The modelled 
Reynolds stress transport equation contains rotation-originated 
redistribution term for the effect of rotation. The modelled 
equation of the turbulence kinetic energy dissipation was modified 
for the effect of streamline curvature. 

4. The effects of the convection and the diffusion were handled 
collectively in the modified Re 3 molds stress model to save the 
computer storage and the computing time. The total effect of the 
convection and the diffusion terms in the Re 3 molds stress transport 
equation was assumed to be related to the production of the 
Reynolds stress. Good agreement between the experimental data 

and the numerical predictions indicates that the above assumption 
is valid for the present problem. 

5. The development of the rotor wake was also predicted with two- 
equation model, which does not include the effect of rotation 
and the streamline curvature. The two-equation model and the 
modified Reynolds stress model show the following difference 

in predictions. The modified Reynolds stress model demonstrates 
substantial Improvements in predicting the mean velocity as well 
as turbulence quantities, compared to the two-equation model. The 
difference in predictions comes mainly from the effect of rotation. 
The two-equation model almost predicts isotropic turbulence 
intensities while the modified Reynolds stress model predicts 
the components of the intensities. Therefore, the two-equation 
model cannot accurately predict turbulent flow whose structure 
is anisotropic. The modified Reynolds stress model predicts 
turbulence Intensity and shear stress components accurately. The 
Increase in radial turbulence Intensity, observed experimentally, 
is also predicted well. The modified Reynolds stress model also 
predicts more accurate shear stresses than the two-equation model. 
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6. Two-dimensional asymmetric turbulent wakes from airfoil were also 
studied numerically in order to Investigate the effect of the 
streamline curvature. Three turbulence closure models were used 
for the prediction. The three turbulence closure models show the 
following differences in predicting two-dimensional asymmetric 
turbulent wakes: 

All three models do not predict the effects of streamline 
curvature accurately when the existing dissipation equation is 
used. At near wake, higher turbulence intensity and shear stress 
are predicted on the pressure side of the wake than are found in 
the experimental results, while lower turbulence intensity and shear 
stress are predicted in the suction side of the wake. 

When the dissipation equation is modified, all three models 
predict mean velocity and shear stress distributions accurately. 

The two-equation model predicts only Isotropic turbulence 
intensities, while the modified Reynolds stress and Reynolds stress 
models predict individual turbulence intensities. 

The modified Reynolds stress and Reynolds stress models 
predict both the mean velocity and turbulence quantities accurately. 
However, computer time and storage can be significantly reduced by 
using the second model. 

7. The effects of rotation and streamline curvature on the development 
of the turbulent wake is substantial and these effects should be 
properly included in the turbulence closure modeling for the 
accurate numerical prediction. 

The prediction scheme reported in this investigation is part 
of a program which will provide a full prediction method for the 
actual flow field inside the turbomachinery. The present prediction 
scheme, combined with the numerical scheme for the three-dimensional 
flow fields inside the rotor blade passages will be invaluable in 
the development of the present status of technology in this field. 
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APPENDIX A 

LISTING OF THE ROTOR WAKE CODE 


C 

C KOTOR WAKK CODi. 

C 

C NUMERICAL PREDICTION OF TURBULENT ROTOR WAKES BY 

0 MODIFIED RSH MODEL WHICH INCLUDES EFFECTS OF 

C STREAMLINE CURVATURE AND ROTATION. 

C 

C 

C INPUT DATA CARD FORMAT AND BRIEF EXPLANATION ARE 

C AS FOLLOWS. 


C 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 


1. FIRST CARD - lAL , DEWS , SOR( 1 5 , 2F 10 . 4 J 

lAL; NUMBER OF ITERATION 
DENS: DENSITY OF FLUID 

SOR; 30R FACTOR 

2. SECOND CARD-DR, RS,DN(3FiO. 3) 

FINITE DIFFERENCES IN R,S,N 
DIRECTIONS 

3. THIRD CARD-EF ,DF , SF(3F10. 3 ) 

EF ; FREE-STREAM TURBULENT 
KINETIC ENERGY 
DF: FREE-STREAM TURBULENT 

ENERGY DISSIPATION RATE 
SF : DUMMY PARAMETER 

4. FOURTH CARD-0M(F10.3 j 

OM: ANGULAR SPEED OF ROTOR 

5 FIFTH CARD-AL(Fi0.3) 

AL: RADIAL VARIATION OF 

FLOW OUTLET ANGLE 

6 SIXTH CARD-VNON( 1) , B( 1 ) , SA( 1 ) ,R( 1) 

(4F10.3) 

VNON(l): FREE-STREAM VELOCITY 
AT HUB 

B(l). OUTLET ANGLE AT HUB 
SAC 1 ) : SPACING AT HUB 

' R(l): R-COODINATE OF HUB 

7 SEVENTH CARD-BC i , 3 AC 1 , VNONC 1 ( 3F i 0 , 3 ) 

BCl: RADIAL VARIATION OF OUTLET 

ANGLE 

SACl: RADIAL VARIATION OF 

SPACING 

VNONCl: RADIAL VARIATION OF 

FREE-STREAM VELOCITY 

8 EIGHTH-51ST CARDS SPECIFY BOUNDAKY CONDITIONS 
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C FOR ELLIPTIC CALCULATION AND DEPEND ON 

C CURVE-FITTING. 

C 

C 

C THE RESULTS ARE PRINTED AND AUTOMATICALLY 

C STORED IN RJE BAT-FILES FOR NEXT ITERATION. 

C 

C 

C 

C 

c 

DIMENSION U(1 i , lOO; ,V(11 , 100) ,W(1 1 , 100) ,UN(7, 100) 
X,N10(10, 100) ,VN(7 , 100) ,WN(7 , 100) ,B(20) ,A(4, 300) 

X, SA(20) , UL(7 , 100) , VLL(7 , 1 00 ) , WL ( 7 , 1 00 ) , VN0N(20) 

X,ENL(7 , 100) ,DIL(7 . 100) ,P(7 , 100) ,PL(7 , 100) ,PN(7 , 100) 
X,N101 (20) ,EN( 1 1 , 100) ,DI( 1 1 , 100) ,G(4 ,300) 

X,ENN(7 , 100) ,DIN(7 , 100 ) , U1 ( 100) ,W1 (100) ,EN1 ( 100) ,DI1(100) 
X,R(20) , IMAX(20) , XC ( 4 , 1 0 , 1 1 ) , RS ( 6 , 1 00 ) , AR ( 6 , 7 ) 

C 


READ(5,27) IAL,DENS,SOR 
2 7 FORMAT Cl 5, 2F1 0.4) 

READ(5,41) DR,DS,DN 

41 FORMA'K 3F 10 . 3 ) 

4991 FORMAT(4F10.3) 

DN0=DN 

READ (5, 41) EF,DF,SF 
C ANGULAR VELOCITY 

READ (5, 42) OM 

42 FORMAT(F10. 3) 

C INVERSE OF ANGLE VARIATION 

READ(3,42) AL 

READ (5 ,4991) VNON(l ) , B ( 1 ) , SA( 1 ) , R ( 1 ) 
READ(3,41) BCl , SAC 1 , VNONC 1 
DO 61 1=2,7 

B(I)=B(1 )+FLOAT( I-l )*BC1 
R(I)=R(1 )+FLOAT( I-l )*DR 
SA( I )=SA( 1 )+FLOAT( I-i )*SAC1 
VNON( I )=VNON( 1 )-FLOAT(I-l )*VN0NC1 
61 CONTINUE 

VNON(d ) = VNON( 1 ) 

VNON(9 ) = VNON( 1 ) 

VNON(10)=VNON(7) 

VNON(l 1 )=VNON(7) 

SA(8)=SA(1) 

SA(9)=SA( 1 ) 

SA(10)=SA(7) 

SA(11)=‘SA(7) 

C NUMBER OF N 

IMAX( 1 )=75 
DO 330 1-1,7 

IMAX(I)-IMAX(1 )+ 2*(I-1) 

330 CONTINUE 



o rj r> o n n o 


142 


IMAX<d)»IMAX(l ) 

IMAX(9)-IMAX(1 ) 

IMAX(10)»IMAX(7) 

IMAX(11)»IMAX(7) 

GETTING INITIAL VELOCITY DISTRIBUTION 

DATA 2 ******** 

DO 301 1-1,11 

READ(5,351)XC(1 , 1 , I ) , XC( 1 , 2 , 1 ) , XC( 1 , 3 , 1 ) , XC ( 1 , 4 , 1 ) ., 
$XC(1,5,I),XC(1 ,6,I),XC(1,7,I) 

301 CONTINUE 

351 FORMAT(7F7 .4 ) 

DO 310 1=1,11 

READ(5 , 352)XC(3 , 1 , 1 ) , XC ( 3 , 2 , X) , XC( 3 , 3 , I ) , XC ( 3 , 4 , I ) , 
$XC(3,5,I),XC(3,6,I),XC(3,7,I),XC(3,a,I),XC(3,9,I) 

310 CONTINUE 

352 FORMAT(9F6.3) 

DO 320 I»l,ll 

READ(5,305)XC(2,1,I),XC(2,2,I),XC(2,3,I),XC(2.4,I) . 
$XC(2,5,I) ,XC(2,6.I),XC(2,7,I),XC(2.8,I) 

320 CONTINUE 
305 fORMAT(8Fb.3) 

DO 307 1-1,11 

READ(5,30U)XC(4, 1 ,1) ,XC(4,2, 1) ,XC(4,3, I) ,XC(4,4,I) , 
$XC(4,5,I),XC(4,6,I),XC(4,7,I),XC(4,8,I),XC(4,9,I). 
$XC(4 , 10, I) 

307 CONTINUE 

308 FORMAT ( 10F6. 3) 

IF(IAL,NE.l) GO TO 3600 
DO 309 1-1,4 

IF( I .£Q, 1 )Jl-7 
IF(I»EQ.2)Jl-8 
IF(I.EQ.3)Jl-9 
IF(I.EQ.4)J1»10 
DO 209 J-l,Jl 
DO 219 K=2,5 

XC(I, J,K)-XC(I, J, 1 )+FL0AT(K-l )*(XC(X , J,6)-XC(I , J, 1) )/ 
$5.0 

219 CONTINUE 
209 CONTINUE 

309 CONTINUE 

INITIAL DISTRIBUTION OF UNKNOWNS 

DO 335 1-1,11 
IMM-IMAX(I) 

IM-(IMM+i)/2 
DO 336 J-1,IMM 
IMMl-IM-20 
IMM2-IM+20 

IF(J .LT. IMMl .OR. J.GT.IMM2)G0 TO 361 
SN»FLOAT(J-IM)*DN/4 .0 
GO TO 365 
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•361 

SN-FL0AT(IMMl-IM)*DN/4 .0 
X+DN*FL0AT( J-IMMl ) 

365 Y1PS“XC(1 ,A,I) 

IFCJ.LT. 1M)Y1PS»XC(1 ,3,1) 

E10»-0.6931*(SN/Y1PS)**2 
IF(K10. LT.-i0.p)E10 — 10.0 
SN3»SN/Y1PS 

U(I , J)-VNON(l)*XC(l , 1 , I)"VNON(I)*XC(l , 2 , I ) *EXP( E 1 0 ) 
X-VNON(I)*XC(l ,2, I)*(“l .0+SNS*XC(l , 5 , I )+XC ( 1 , 6 , I ) *SNS**2 
$+XC( 1 , 7 * I)*SNS**3) 

Y3 = oN/}>A(I)*2 .0 
IF(J.LT.IM) GO TO 215 
IF(Y3.GE.XC(3,2,I)> GO TO 220 
VL=VN0N(I)*(XC(3 , 3, I)-XC(3, 1 , I) )/XC(3,2, I) 

W(I , J) = VL*Y3 + VNON(I)*XC(3 , 3 , I)-VL*XC(3 , 2,1) 

GO TO 325 

220 IF(Y3.GE.XC(3,4,1)) GO TO 225 

VE=VN0N(I)*(XC(3 > 5 , I ) -XC( 3 , J , I ) ) / ( XC ( 3 , 4 , I ) -XC ( 3 , 2 , 1 ) ) 
W( I , J)=VL*Y3+VN0N( I) v>t.XGC3 , 5 , I)-VL*XC( 3 ,4,1) 

GO TO 325 

225 W(I , J)=VN0N(I)*XC(3, 5, I) 

GO TO 325 

2,15 IF( Y3 . LE .XC(3 , 6 , I ) ) GO TO 226 

VL=(XC(3 , 7 , I)-XC(3 , 1,1) )*VN0N(1)/XC(3 . 6,1) 

W(I,J)“VL*Y3 + VNON(I)*XC(3,7 . I)-VL-^XC(3,6, I) 

GO TO 325 

226 IF( Y3 . LE .XC(3 , o , 1) ) GO TO 229 

VL=<XC(3,9, I)-XC(3,7 , 1) ) * VNON ( I ) / C XL ( 3 , 8 , I)-XC(3,6, I)) 
W( 1 , J) = VL*Y3 + VNON(i)*XCC J , 9 , 1)-VL*XC(3 , 8,1) 

GO TO 325 

229 W(I , J)=VNON(I)*XC(3 , 9 , 1) 

325 IF( J.LT. Il'l) GO TO 240 

IF(Y3 .Gf .XC(2 , 1 , I) ) GO TO. 24 5 

V(I , J)=XC(2 ,2 , I)/XC(2 , 1 , I)*Y3*VN0N(I) 

GO TO 336 

245 IF (Y3.GT.XC(2,3, I) ) GO TO 247 

VL = VNON(I)*(XC(2,4,l)-XC(2,2,I))/(XC(2,3,I)-XC(2,l ,D) 
V(I,J)«=VL*Y3+XC(2,4 , I ) *VNON ( I ) - VL*XC ( 2 , 3 , I ) 

GO TO 336 

247 V(I,J)=VNON(I)*XC(2,4,I) 

GO TO 336 

240 XF(Y3.LE. XC ( 2 , 5 , I ) ) GO TO 250 

V(I , J)=VNON(I)*XC(2 . 6 , I)/XC(2 , 5, I)*Y3 
GO TO 336 

250 IF( Y3 . LE .XG(2 , 7 , 1) ) GO TO 255 

VL = VNON(I)*(XC(2,8,I)-XC(2,'6,I))/(XC(2,7,I)-XC(2,5,l)) 
V(1 , J)-VL*Y3+VNON(I)*XC(2 ,a,I)-VL*XC(2, 7 ,1) 

GO TO 336 

235 V(I,J)«VNON(I)*XC(2,8,I) 

336 CONTINUE 

IUC-1MM-50 
IMCl -IMC+1 
DO 410 J3«l,50 
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U(I,J3)«U(I,51) 

410 CONTINUE 

DO 415 J3-IMC1.IMM 
U(I , J3)=U(I , IMC) 

415 CONTINUE 
335 CONTINUE 

DO 3335 I«l,li 
IMM=IHAX( I ) 

IM=( IMM-Hl )/2 

IHMl=IM-20 

1MM2-IM+20 

DO 3336 J=»i . IKM 

IMX=lMHl-30 

lMY = IHH2-r3U 

lECJ.LT. IMX.OR. J.GT . IMY) GO TO 5oO 
1K(J.LT.IMM1 .0R.J.GT.IMM2) CO TO 3361 
DN=DN0/4 . 0 
SN=r LOAT( J-IM) *0N 
GO TO 3365 

3361 IF(J.GT.IM) Iimi=lilM2 
DN^DNO 

S N = F LOAT ( IMM 1 - IM ) *DN / 4 . 0+DN f LOAT C J - 1 lill i ) 

3365 Y3=SN/3A(I)*2 , 0 

IF(J.EQ. 1 .OK.J.FQ. IMM) GO TO 5<iu 
ECE=0 .007 

IF(J.LE.IH) cCE=0.002 
IF(I.EQ.7) ECE=ECE*5,U 

TT = -XC(4 , 2 , l)*VNON(I )**2/(XC(4 , 1 , I) *'i . 2 ) 

UG=(U(I , J+1 )-U( I , J-1 ) )/2 .0/DN 
WG=(W(I , J+1 )-W(I , J-I ) )/2.0/DN 

IF(J.EQ, IMMl)UC = 0.5*((U(l , J-rl ) -U (1 » J ) ) / DN-r ( U ( I ,J) 
))/4.0/DN) 

IF ( J . eg . IMMl ) WG = 0 . 5* ( ( W( I , J + i )-W( I , J ) > / DW-f ( W( I , J ) 
$-W(I,J-l ))/4.0/DN) 

IF( J ,EQ. IMH2)UG = 0 . 5 ■’« ( ( U ( I .J+1 ) -U ( i , J ) ) /4 . 0 /D W+ ( U ( i , J ) 
$-UCi. J-1 ) )/DN) 

IF( J . EQ. IMM2)WC = 0. j-K (W( I , J-r 1 ) - W ( I , J ) ) / 4 . 0/DN+(W( I , J ) 
$-W(I , J-1 ) )/DN) 

EN(I, J) = ECE*(ABS(UG)-t-0. 5»ABS(VJG) )/3.0 
DIC = O. 60-0.075*FLOAT(I-i ) 

DIC=DIC*0 . 0 

IF(I .EQ. 6 )DIC=U . 3 

Diri . J)=l •0*DIG*‘EN(1, J)*AB3(UG) 

$+0.3*DIC*LN(I , J)*ABS(WC) 

IF(EN(I,J) .GT.EF .AND.DI(I,J) .GT.DF) GO TO 3336 
530 DI(I,J)=DF 
LN(I , J)=EF 
3330 CONTINUE 

AAA AAA A vV A 

IF(I.NE.2) GO TO 3335 

WRITE(6,691) (J,EN(I , J) ,DI(I , J) , J=i , IMM) 

691 FORMAT(4(' ',13,' ',2F10.1)) 

3335 CONTINUE 
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C 

C 

C 

C 

0 hlt-k *** 

C 

C 

DO 451 1-1,7 
IMM-IMAX(I) 

IM-(IMM+1 )/2 
IMMl-IM-20 
IMM2-IM+20 
DO 455 J-2,100 
N10(I , 1 )-l 

N10(I, J)-N10(I,J“1)+10 
IF(N10(I,J ).GT.IMM1) GO TO 461 
455 CONTINUE 

461 N10(I , J)-IMM1 

462 J-J+1 

N10(I,J)»N10(I,J-1 )+2 
IF(N10(I , J) .GT . IMM2)GO TO 463 
GO TO 462 

463 N10(I , J)=N10(I , J-1 )+10 

464 J=J+1 

N10(I, J)»N10(I,J-1 )+10 
IF(N10(I , J) .GT. IMM)GO TO 466 
GO TO 464 
466 N101(I)=J-1 
N20-N101(I) 

WRITE(6, 1421 ) I 

1421 F0RMAT(1H , 'INITIAL CONDITION AT STATION- 13 ) 

WRITE (6, 4 70) (N10(I,J),U(I,N10(I.J)).V(I,N10(I,J) 
X) ,W(I ,N10(I ,J)),EN(I,N10(I, J) ) ,DI(I ,N10(I, J>), 
$J«1 ,N20) 

470 FORMAT(I5,5F10.2,5H ,15,5FI0.2) 

451 CONTINUE 

DO 1820 1-1,7 
KM-IMAX(I) 

DO 2011 J-1,KM 
P(l, J)-~0.5*0M**2*R(I)**2 
2011 CONTINUE 

WRITE (14, 830) (U(I,J),V(I,J),W(I,J),P(I,J),EN(I,J 
$),DI(I,J),J-1,KM) 

1820 CONTINUE 
C 

0 A4rA it it it *** 

C 

C 

C 

DO 700 K-1 , 12 

GET BOUNDARY CONDITION FOR R-SWEEP 
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IMM-IMAX(l) 

FK-0.0 

IF(K.EQ.l) GO TO 277 

fk-cfloatU“1)) 

277 DO 1320 

U(1 , J)-(U(9, J)-U(8, J))/ll .0*FK+U(8, J) 

V( 1 , J)-(V(9, J)-V(8, J))/l 1 .0*FK+V(8, J) 

W( 1 , J)-(W(9, J)-W(8, J))/ll .0*FK+W(8, J) 

EN(1 . J)-(EN(9, J)-EN(8, J) )/l 1 .0*FK+EN(8, J) 
Did . J)-(DI(9, J)-DI(8, J))/l 1 .0*FK+DI(8, J) 
1320 CONTINUE 

IMM7-IMAX(7) 

DO 71320 J-1 ,IMM7 

U(7,J)-( U(11,J)- U(10,J))/11.0*FK+U(i0,J) 
V(7,J)-( V(11,J)- V(10, J) )/li •0*FK+V(i0. J) 
W(7,J)"( W(11,J)“ W(10,J))/11.0*FK+W(10,J) 
EN(7 , J)-(t:N( 1 1 , J)-EN(iO, J))/ll .0*FK+£N(10, J) 
DI(7 , J)=(D1(1 1 , J)-Dl(lO.J) )/11.0*FK+DI( 10, J) 
71320 CONTINUE 


STARTING ITERATION 


SS»DS*FLOAT(K) 

DO 10 1-2,6 
IMM-IMAX(I) 

IMMOl-IMM-1 
IM-(IMM+l)/2 
IMML»IMAX(I-1 ) 

IML-CIMML+I )/2 
IMMl-IM-20 
IMM2-IM+20 
S-SIN(B(D) 

C-COS(B(I) ) 

ID-IM-IML 

CH1»-SS*S**4/R(I)**2 

CH2--S**2/R(I) 

CH3-SS**2/2.0/R(I)**3*2 ,0*S**4-SS**2/R(I)**2/AL*S**3*C 
CH4-SS**2/R(I)**3*S*C-1 . 0 /AL-S S**2/R( I ) **2 / AL* ( S*C ) **2 
CHS — SS*S**3*C/R(I)**2 
CH6-CH5+SS/R(I)/AL*S**2 

CH7-SS**2/R(I)**3*S**3*C+1 .0/AL-SS**2/R(I)**2/AL 
X*s**2*(C**2+l .0) 

CH8--S*C/R(I) 

CH9-SS/R(I)**2*S*C-SS/R(I)/AL*C**2 

CH10-SS/R(I)**2*S**2-SS/R(I)*S*C/AL 

CH11»-*SS**3/R(I)**4*S**4-2.0*SS**3/R(I)**3/AL*S**3*C 
CH12— SS**2/R(I)**3*S**2+2,0*SS**2/R(I)**2/AL*S*C 
CH14--SS**3/R(I)**4*S**3*C+SS**3/R(I)**3 
$/AL*S**2*(i .0+2.0*C**2) 

CH15— SS/R(I)**2*(S*C)**2+SS/R(I)/AL*S*C 
CH16-SS**2/R(I)**3*S**2*C**2-SS**2/R(I>**2/AL*S*C 
$*(C**2+1 .0) 

AA— SS*S/R(I)+SS*C/AL 
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DO 3 

U1(J )-U(I,J) 

W1(J)-W(I»J) 

EN1(J)»EN(I,J) 

DI1(J)-DI(I, J) 

3 CONTINUE 

DO 905 J-2.IMM01 
DO 910 N-1,6 
DO 915 M-1 ,6 
AR(N,M)-0 .0 
915 CONTINUE 

AR(N,7)-0.0 
910 CONTINUE 

DN-.DN0/4 .0 

IF( J ,LT. IMMl .OR. J.GT.IMM2) DN-ONO 
GUS-0.0 

GU-(U(I , J+1 )-U(X , J-1 ) )/2.0/DN 
GW-(W(I , J+1 )-W(I , J-i ) )/2.0/DN 
TE»EN(I , J) 

TD=DI(I , J) 

AR( 1 , 1 )»-l .5*TD/TE-4 .0*OM*S**3*3S/R(I)'-GUS-U(1 , J)*CH1 
X-W( I , J)*CH3 

AR( 1 ,4)— 0 .8*(-SS/R(I)*S*C+SS/AL)*(GU+U(I . J)*CU5+ 

$W(I » J)*CH4 

X)+4 ,0*0M*S*( 1 .0+SS**2/R(I )**2*S**2-2 .0*SS**2/R(I) 
$-^*2/AL*S*C) 

AR( 1 ,5)--0.8*(GU+U(I , J)*CH1+W( I , J)*CH3)-4.0*OM*S 
$**2*C*SS/R(I) 

AR( 1 , 7 )— TD/3 .0*2 . 2*( 1 .0+SS**2/R(I)**2*S**4 ) 

AR(2 ,2)— 1 .5*TD/TE+4 .0*OM*C*(-SS/R(I)*S*C+SS/AL) 

AR(2 ,5)»4 .0*OM*C*(-SS/R(I)*3**2 ) 

AR(2 , 6)»4 .0*OM*C*(1 . O+S S* * 2 / R ( I ) * * 2 * S * *2- 2 . 0 *S S * * 2 
$/R(I)/AL*S*C) 

AR(2 , 7 )— TD/3 .0*2 . 2*( 1 . 0+SS**2/R(I ) **2*S**2*C**2 
X-2,0*SS**2/R(I)/AL*S*C) 

AR(3 , 3)“-0 .8*(-SS/R(I)*S*C+SS/AL)*(GW+U(I , J)*C118 
$+W(I , J)*CH9) 

$-0.8*(U(I , J)*CH2+W(I ,J)*CH10)*SS/R(I)*S**2+4 .0*OM 
$*S**3*SS/R(I) 

AR(3 , 4 )=-4 .0*OM*S-0 . 8* (U( I , J ) *CH2+W( I , J)*CH10) 

AR(3 , 6)»-0 .8*(GW+CH8*U(I, J)+W(I , J)*CH9)-4.0*0M*C 
AR(3»7) — TD*2. 2/3.0 

AR(4 , 1 )— 0 .4*(U(I , J)*CH2+W(I , J)*CH10)-2.0*OM*S 
AR(4 , 3)— 0 .4*(-SS/R(l)*S*C+SS/AL)*(GU+U(I , J)*CH5 
$+W(I , J)*CH4) 

$+0.4*SS/R(I)*S**2*(GUS+U(I , J)*CH1+W( I , J)*CH3)+2 .0 
$*OM*S*( 1 .0 

$+SS**2/R(I)**2^S**2-2.0*3S**2/R(l)/AL*S*C) 

AR(4 »4) — 0.4*(GUS+U(I , J)*CH1+W(I , J)*CH3)“0 .4*(-SS/ 
$R(I)*S*C 

§+SS/AL)+0.4*SS/R(I)*S**2*(U(I, J)*CH2+W( I , J)*CH10) 
$“2.0*0M*C*(-SS/R(I)*S*C+SS/AL)-1 .5*TD/TE 
AR(4,5)— 0.4*(GW+U(I , J)*CH8+W(I, J)*CH9)-2.0*0M*C 
AR(4,6>— 0.4*(GU+U(I , J)*CH5+W(I »J)*CH4)+2.0*0M*S 
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$*(-SS/R(I)*S*C 

$+SS/AL) 

AR(4 , 7 )«=2.2/3.0*TD*SS/R(I)*S**2 
AR(5 , 1 )=-2 .0*OM*C*SS/R(I)*S**2 

AR(5,2)=-0.4*(GU+U(I ,J)*CH5+W(I, J)*CH4)+2.0*OM*S* 
$(-SS/R(I)*S*C 
$+S S/ AL ) 

AR(5 ,4)=2 .0*OM*C*(1 .0+SS**2/R(l)**2*S**2) 

AR(5 , 5)=-0.4*(GUS+U(I , J)*CH1+WU . J)*CH3)-2.0*OM 
$*S**3*SS/R(I) 

$+2 .0*0M*C*(-SS/R(I)*S*C+SS/AL)-1 . 5*TD/TE 
AR(5 , 6)=-0.4*(-SS/R(I)*S*C+SS/AL)*(GU+U(I , J)*CH5 
$+W(I , J)*CH4) 

$+2 .0*0M*S*(1 .0+SS**2/R(I)**2*S**2-2 .0*SS**2/R(I)/ 
$AL*S*C) 

AR(5 , 7)=-TD*2 .2/3.0*(SS**2/R(I)**2*C*S**3-SS**2/ 
$R(I)/AL*S**2) 

AR(6 , 2)=-0.4*(GW+U(I , J)*CH8+W(I , J)*CH9)-2.0*0M*C 
AR(6,3)=2 .0*0M*C*(1 . 0+SS**2/R( I ) **2*S**2 
$-2 .0*SS**2/R(I)/AL*S*C) 

AR(6 ,4)=-2 .0*0M*C*SS/R(I)*S**2 

AR(6 ,5)=-0.4*(U(I , J)*CH2+W(I, J)*CH10)-2.0*OM*S 
AR(6 , 6)=-0 ,4*(-SS/R(I)*S*C+SS/AL)*(GW+U(I , J)*CH8 
$+W( I , J)*CH9) 

$+0 ,4*SS/R(I)*S**2*(U(I , J)*CH2+W( I , J)*CH10)~2 .0*0M 
$*C*(-SS/R(I) 

$*S*C+SS/AL)+2 ,0*0M*S*SS/R(I)*S**2+2 .0*0M*C*(-SS/R 
$ (I)*S*C+SS/AL) 

$-1 .5*TD/TE 

AR(6 , 7 )*2 .2/3.0*TD*SS/R(I)*S**2 

DO 930 II»1,5 

IJ=II+1 

AR(II,II)=1 .0/AR(II,II) 

AR(II, 7 )=AR(II, 7)*AR(II, II) 

DO 931 N=IJ,6 

AR(II,N)=AR(II,N)*AR(II, II) 

931 CONTINUE 

DO 932 N»IJ,6 

AR(N,7)=AR(N,7)-AR(II, 7)*AR(N, II) 

932 CONTINUE 

DO 933 N=IJ,6 
DO 934 M=IJ,6 

AR(M,N)=AR(M,N)-AR(M, II)*AR(II,N) 

934 CONTINUE 

933 CONTINUE 
930 CONTINUE 

AR(6 , 7 )=AR(6 , 7)/AR(6 , 6) 

NV«6 

943 NV-NV-1 

IF(NV.EQ.O) GO TO 944 

NX=NV+1 

DO 936 N=NX,6 

AR(NV,7)-AR(NV,7)-AR(NV,N)*AR(N,7) 

936 CONTINUE 
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60 TO 943 

944 RS(1,J)-AR(L,7) 

RS(2,J)-AR(2,7) 

RS(3,J)-AR(3,7) 

RS(4,J)-AR(4,7) 

RS(5, J)-AR(5,7) 

RS(6, J)«AR(6,7) 

PR^— (GU*RS(5, J)+6W*RS(6, J))-TD 

WR1TE(6,6330) 1 , J , 6U , 6W ,RS ( 1 , J) ,RS (2 » J) ,RS (3 , J) ,RS 
$(4 , J) ,RS(5, J),RS(6,J) ,PRE 

6330 FORMAT(' ',214,' ',2F10.2,' ',6F10.4,' ',F12.1) 

90S CONTINUE 
MM-0 

15 MM-MM+1 
NEW S 

PR-2.0*0M*U(I,IMM)*S+R(I)*0M**2-U(I,IMM)**2*CH2 
PR1-U(I , IMM)**2*CHl+2 .0*0M*(-SS*S/R(I)+SS*C/AL)*S**2* 
$U(1,IMM) 

DO 20 J-1 , IMM 
RUV-0.0 
IMMl 1-lMMl+l 
IMM22-IMM2-1 

IFCJ.LE.lMMll .OR. J.6E.1MM22) GO TO 965 
DN-DNO/4 .0 

RUV-(RS(5, J+1)-RS(5 , J-1 ) )/2 .0/DN 
965 IF((J-ID).LT.l) GO TO 145 

IF( (J-ID) .GT.IMML) GO TO 146 
UNN-UCI-I , J-ID) 

60 TO 141 

145 UNN-U(I-1,1) 

GO TO 141 

146 UNN-U(I-1 , IMML) 

141 DN«DNO 

IF( J.GE.lHMl .AND. J.LE.IMH2)DN-DN0/4 .0 
A(1 , J)— V(I, J)/2.0/DN 
A(2 , J)-U(I , J)/DS+W(I , J)/DR 
A(3, J)-V(l,J)/2.0/DN 

A(4 , J)-PR1-U(I , J)**2*CHl+2 .0*OM*(-SS*S/R(I)+SS*C/AL 
X)*U(I . J)*S**2"RUV+U(I, J)*U1(J)/DS 
X+W(I, J)*UNN/DR 
1F(J.EQ.1MH1)60 TO 19 
1F(J.£Q.IMM2)G0 TO 18 
GO TO 20 

19 A(1 , J)— V<I,J)/8.0/DN 

A(2 , J)-U(I, J)/DS+W(I , J)/DR-3 .0/8.0*V(I, J)/DN 
A(3, J)-V(I, J)/2.0/DN 
GO TO 20 

18 A(1 , J)— V(I, J)/2.0/DN 
A(3,J)-V(I,J)/8.0/DN 

A(2 ,J)»U(I,J)/DS+3.0/8.0*V(I, J)/DN+W(I ,J)/DR 

20 CONTINUE 

A(4,1)-A(4,1)-A(1,1)*U(I,1) 


( 
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A(4 , 1MM)-A(4 , 1MM)-A(3 , IHH)*U( I , IMH) 
CALL ELIM(A»IMM) 

DO 101 J"1,IMH 
UN(I, J)-A(4, J) 

101 CONTINUE 


NEW R 

DO 30 J-1 »1HM 
RWV-0.0 
IMMll-lMHl+1 
IMM22-IMM2-1 

1F(J.LE.1MM11.0R.J.GE.IMH22) GO TO 966 
DN-DNO/4.0 

RWV»(RS(6, J+1 )-RS(6, J-1) )/2 .0/DN 
966 IF( (J-ID) .LT. 1) GO TO 165 
IF( (J-ID) .GT.IMMDGO TO 166 
WNN- W(I-1 , J-ID) 

GO TO 161 

165 WNN- W(I-1 , 1 ) 

GO TO 161 

166 WNN- W(I-1 , IMML) 

.161 DN-DNO 

IFCJ.GE .IMMl .AND. J.LE.IMM2)DN-DN0/4.0 
A(1 , J) — V(I, J)/2.0/DN 
A(2, J)-UN(I, J)/DS+W(I, J)/DR 
A(3, J)-V(I,J)/2.0/DN 

A(4 , J)-UN(I, J)*W1 ( J)/DS+W(1, J)*WNN/DR-UN(I , J) 
X**2*CH2-RWV-2.0*OM*UN(I , J)*S+R(I)*0M**2-PR 
IF(J.EQ.IMM1)G0 TO 29 
1F(J.£Q.1MH2)G0 TO 28 
GO TO 30 

29 A( 1 , J)— V(I , J)/8.0/DN 
A(3, J)-V(I,J)/2.0/DN 

A(2 »J)-UN(I, J)/DS+W(I,J)/DR-3.0/8.0*V(I, J)/DN 
GO TO 30 

28 A(1 ,J) — V(I,J)/2.0/DN 
A(3 . J)-V(I , J)/8.0/DN 

A(2 , J)-UN(I , J)/DS+3.0/8.0*V(I » J)/DN+W(I , J)/DR 

30 CONTINUE 

A(4,1)-A(4,1)-A(1,1)*W(I,1) 

A(4, IMM)-A(4,IMH)-A(3, IHM)*W(I, IMM) 

CALL ELIM(A,IHM) 

DO 102 J«1,IMM 
WN(I, J)-A(4, J) 

102 CONTINUE 
UMl-200.0 
DO 103 J-1, IMH 
U(I,J)-UN(I,J) 

W(I,J)-WN(I,J) 

IF(J.LT.IMH1 .OR. J.GT.IMH2)GO TO 103 
IF(U(I , J) .GE .UMl >60 TO 103 
UM1-U(I,J) 

NT-J 
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103 CONTINUE 

1F(MH.GE.2)G0 TO 600 
GO TO 15 

NEW N 

600 NTl-NT+1 

VN(i,NT)-0.0 
DO 40 J«NT1,1HH 
DN-DNO 

lF(J.LE.IMH2)DN«DN0/4.0 
IF( (J-ID) .GT.IMMDGO TO 176 
WNN- W(I-1 , J-ID) 

GO TO 171 

176 WNN- W(I-1 . IMML) 

171 GM— ( U(l, J)-Ui(J))/DS 
VN(I , J)-VN(I , J-1 )+GM*DN 
40 CONTINUE 
NT2-NT-1 
DO 640 J-1 ,NT2 
DN-DNO 
JJ-NT-J 

IF( JJ.GE. IMMl ) DN-DNO/ 4.0 
IF( ( JJ-ID) .LT. 1 ) GO TO 175 
WNN-W(I-1 , JJ-ID) 

GO TO 671 
175 WNN-W(I-l.l) 

671 GM— (U(I , JJ)-U1 (JJ))/DS 

VN(I , JJ)-VN(I , JJ+1 )-GM*DN 
640 CONTINUE 

DO 46 J-1 , IMM 
V(I, J)-VN(I,J) 

46 CONTINUE 

DO 1020 J-1 , IMM 
IF((J-ID) .LT.l) GO TO 1145 
IF( (J-ID) .GT. IMML) GO TO 1146 
ENNN- EN(I-1,J-ID) 

DINN- DI(I-1,J-ID) 

GO TO 1141 

1145 ENNN-EN(I-1 , 1 ) 

DINN- DI(t-l,l) 

GO TO 1141 

1146 £NNN-£N(I-1 , IMML) 

DINN- DI(I-1 , IMML) 

1141 DN-DNO 

IF(J.G£.IMM1.AND. J.LE.IMM2) DN-DNO/4.0 

IF(J.EQ.l) GO TO 73 

IF(J.EQ.IMM) GO TO 74 

IF(J.EQ.IMMi) GO TO 75 

IF( J.EQ. IMM2) GO TO 76 

Q1-(U(I,J+1 )-U(I ,J-l))/2.0/DN 

Q2-(W(I , J+1 )-W(I , J-1 ))/2.0/DN 

GO TO 77 

73 Q1-(U(I,J+1)-U(I,J))/DN 
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Q2-(W(I,J+1)-W(I,J))/DN 
GO TO 77 

74 Q1-(U(I,J)-U(I,J-1))/DN 
Q2-(W(I,J)-W(I,J-1))/DN 
GO TO 77 

75 Q1»0.5/DN*(U(I,J+1)-0.75*U(I, J)-0.25*U(I.J-i)) 
Q2»0.5/DN*(W(I.J+1 )-0.75*W(I, J)-0.25*W(I,J-1 ) ) 
GO TO 77 

76 Ql-0.5/DN*(0.25*U(I , J+1 )+0 . 75*U(I , J)-U(I » J-1 ) ) 
Q2-0.5/DN*(0.25*W(I, J+1 )+0 . 7 5*W( I , J ) -W( I , J- 1 ) ) 

77 E-0.09*EN(I,J)**2/DI(I, J) 
PR=E*(Q1**2+Q2**2+Q1*(CH5+CH6)+Q2*(CH7+CH8)) 

A(1 , J)-=-VN(I, J)/2.0/DN-E/DN**2 

A(2 , J)“UN(I , J)/DS+WN(I , J)/DR+E*2 ,0/DN**2 
A(3 , J)=VN(I , J)/2 .0/DN-E/DN**2 
A(4 , J)=UN(I , J)*EN1 ( J)/DS+WN(I , J)*ENNN/DR 
$+PR 

$-Dl(I,J) 

G( 1 , J)=-VN(I , J)/2 .0/DN-E/l .3/DN**2 
G(2 , J)=UN(I , J)/DS+WN(I , J)/DR+E*2 .0/1 .3/DN**2 
G(3 , J)=VN(I , J)/2 .0/DN-E/l .3/DN**2 
G(4, J)=UN(I , J)*DI1(J)/DS+WN(I, J)*DINN/DR 
$ + 1 .44*DI(I,J)/EN(I ,J)*PR 
$-1.92*DI(I,J)**2/EN(I,J) 

IF(J.EQ.IMMl) GO TO 1019 
IF(J.EQ. IMM2) GO TO 1018 
GO TO 1020 

1019 A(i , J)=*-VN( I »J)/8. 0/DN-E/l 0.0/DN**2 

A(2 , J)=UN(I,J)/DS-3.0/8.0*VN(I , J)/DN+WN(I , J)/DR 
$+0 . 5*E/DN**2 

A(3, J)=VN(I,J)/2.0/DN-E/2.5/DN**2 
A(4 , J)=UN(I, J)*EN1 (J)/DS+WN(I » J)*ENNN/DR 
$+PR -DI(I,J) 

G(1 ,J)=-VN( I ,J)/8 .0/DN-E/l 3. 0/DN**2 
G(2 , J)=UN(I , J)/DS-3 .0/3.0*VN(I , J)/DN 
$+WN(I , J)/DR+1 .0/2.6*E/DN**2 
G(3 , J)=VN(I , J)/2 .O/DN-E/2 . 5/1 . 3/DN**2 
G(4, J)=UN(I »J)*DI1 (J)/DS+WN(I, J)*DINN/DR 
$+1 .44*DI(I,J)/EN(I,J)*PR 
$ -1 . 92*DI(I , J)**2/£N(I , J) 

GO TO 1020 

1018 A(1 , J)=-VN(I,J)/2.0/DN-E/2.5/DN**2 
A(3, J)=VN(I , J)/8.0/DN-E/10.0/DN**2 
A(2, J)=UN(I,J)/DS+3.0/8.0*VN(I , J)/DN+WN(I , J)/DR 
$+0 .5*E/DN**2 

A(4, J)=UN(I,J)*EN1(J)/DS+WN(I,J)*ENNN/DR 
$+PR 

$-DI(I,J) 

G(1 »J)=-VN(I,J)/2.0/DN-E/2.5/l .3/DN**2 
G(3 , J)=VN(I , J)/8.0/DN-E/13.0/DN**2 
G(2,J)=UN(I,J)/DS+3.0/8.0*VN(I, J)/DN+WN(I.J)/DR 
$+0.5/1 .3*E/DN**2 

G(4, J)-UN(I, J)*DI1(J)/DS+WN(I » J)*DINN/DR 
$+1 .44*DI(I , J)/EN(I, J)*PR 
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$-1 .92*DI(I, J)**2/liN(I»J) 

1020 CONTINUE 

A(A,1)-A(4,1)-A(1,1)*EN(I,1) 

A(4 , IMM)-A(4 ,IMM)“A(3, IMM)*EN(I,IMM) 

G(4,1)-G(4,1)-G(1,1)*DI(I,1) 

G(4,IMM)-G(4,IMM)-G(3,1MM)*DI(I,IMM) 

CALL EL1M(A,IMM) 

DO 1101 J»1,IMM 
ENN(I , J)«A(4 , J) 

1101 CONTINUE 

CALL EL1H(G,1HM) 

DO 1107 J-1,IMM 
DIN(I.J)»G(4.J) 

1107 CONTINUE 

DO 106 J-1 , IMM 
EN(I , J)=ENN(I , J) 

Did , J)=»DIN(I » J) 

106 CONTINUE 
C 
C 

*** *** *** *** *** 

C 

c 

N20=N10 1 (I) 

WRITE(6 , 2370)K, I 

WRITE(6 ,2470) (N10(I,J),U(I,N10(I,J)),V(I, 

$N10(I, J) ) ,W(I,N10 

X(I,J)),EN(I,Ni0(I,J)),DI(I,N10(I,J)),J-l ,N20) 

2470 FORMAT( I5,5F10.2,5H ,I5,5F10.2) 

2370 FORMATCIH ,'***',' S-S TAT ION* 13 , ' R-STATION 

6=' 13 

10 CONTINUE 

DO 820 1=1,7 
KM=IMAX( I) 

DO 5011 J»1,KM 
F(I, J)=“0.5*OM**2*R(I)**2 
5011 CONTINUE 

WRITE (14, 830) (U(I,J),V(I,J),W(I,J),F(I,J), 

$EN(I, J) ,DI(I, J) , J*1 ,KM) 

820 CONTINUE 

830 FORMAT (3F5. 1 ,F7 . 1 , 2F9 . 1 , 3F5 . 1 . F7 . 1 , 2F9 . 1 ) 

700 CONTINUE 

GO TO 44444 
3600 DO 4700 K*1 , 1 2 

IF(K.NE.l) GO TO 1840 
DO 1830 1=1,7 
KM=IMAX(I) 

READ(5 ,830) (UL( I , J ) , VLL( I , J ) , WL ( I , J ) , PL ( I , J ) , ENL( I , J ) 
$,DIL(I ,J) , J*1 ,KM) 

1830 CONTINUE 

DO 1835 1=1,7 
KM»IMAX(I) 

READ(5,830) (U(I,J),V(I,J),W(I,J),P(I,J),EN(I,J) 

$DI(I , J) , J-1 ,KM) 
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1335 CONTINUE 
1840 DO 1845 1-1,7 
KM-IMAX(I) 

READ(5 ,330)(UN(I ,J) , VN(I , J) ,WN(I , J) ,PN(I, J) ,ENN(1 , J) 
$,DIN(I , J) , J-1 ,KM) 

1845 CONTINUE 


SS-DS*FLOAT(K) 

DO 4010 1-2,6 
IMM-lHAX(l) 

IMMOl-IMM-1 
IM-(IMM+l)/2 
IMML-IMAX( I-l) 

IML-(IMHL+1 )/2 
IMMl-IM-20 
IMM2-IM+20 
S-SIN(B(I) ) 

C-COS(B(l) ) 

ID-IM-IML 

CHI— SS*S**4/R(I)**2 
CH2 — S**2/R(I) 

CH3-SS**2/2 .0/R(I)**3*2.0*S**4-SS**2/R(I)**2/AL*S**3*C 
CH4-SS**2/R(I)**3*S*C-1 .0/AL-SS**2/R(I)**2/AL*<S*C)**2 
CH5--SS*S**3*C/R(I)**2 
CH6-CH5+SS/R(I)/AL*S**2 

CH7-SS**2/R(I)**3*S**3*C+1 . 0/AL-SS**2/R( 1) **2/AL 
X*s**2*(C**2+l .0) 

CH8— S*C/R(I) 

CU9-SS/R(1)**2*S*C-SS/R(I)/AL*C**2 

CH10-SS/R(I)**2*S**2-SS/R(I)*S*C/AL 

CHll— SS**3/R(I)**4*S**4-2 .0*SS**3/R(I)**3/AL*S**3*C 
CH12— SS**2/R(I)**3*S**2+2.0*SS**2/R(I)**2/AL*S*C 
CH14— SS**3/R(I)**4*S**3*C+SS**3/R(I)**3 
$/AL*S**2*(l .0+2.0*C**2) 

CH15— SS/R(I)**2*(S*C)**2+SS/R(I)/AL*S*C 
CH16-SS* ' 2/R(I)**3>‘S**2*C**2-SS**2/R(I)**2/AL*S*C 
$*(C**2+1 .0) 

DO 4905 J-2,1HM01 
DO 4910 N-1,6 
DO 4915 M-1,6 
AR(N,M)-0.0 
4915 CONTINUE 

AR(N,7)-0.0 
4910 CONTINUE 

DN-DNO/4 .0 

IFCJ.LT.IHMI .OR. J.GT.1MM2) DN-DNO 
GUS-(UN(I,J)-UL(I,J))/2.0/DS 
GWS-(WN(I , J)-HL(I , J) )/2.0/DS 
GU-(U(I, J+1 )-U(I,J-l))/2.0/DN 
GW-(W(I , J+1 )-W(I , J-1 ) )/2.0/DN 
TE-EN(I , J) 

TD-DI(I,J) 

AR(1 , 1)— 1 .5*TD/TE-4.0*0M*S**3*SS/R(I)-GUS-U(I, J)*CH1 
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X-W(I,J)*CH3 

AR(1 ,4)— 0.8*(-SS/R(I)*S*C+SS/AL)*(GU+U(I , J)*CH5+ 
$W(I , J)*CH4 

X)+4 .0*OM*S*( 1 ,0+SS**2/R(I)**2*S**2-2 .0*SS**2/R(I) 
$**2/AL*S*C) 

AR( i , 5)— 0 ,8*(GU+U(I , J)*CH1+W(I , J)*CH3)-4 ,0*OM*S 
$**2*C*SS/R(I) 

AR(1 .7)— TD/3.0*2.2*(1 .0+SS**2/R(I)**2*S**4) 

AR(2 ,2)— 1 .5*TD/TE+4 ,0*OM*C*(-SS/R(I)*S*C+SS/AL) 
AR(2 ,5)-4.0*0M*C*(-SS/R(I)*S**2) 

AR(2 ,6)-4.0*OM*C*(l . 0+S S**2 / R( I) **2 *S**2-2 . 0 *S S**2 
$/R(I)/AL*S*C) 

AR(2 ,7 )— TD/3 .0*2 .2*C1 . 0+S S**2 /R( I) **2 *S**2 *C* *2 
X-2 .0*SS**2/R(I)/AL*S*C) 

AR(3 , 3)— 0.8*(-SS/R(I)*S*C+SS/AL)*(GW+U(I , J)*CH8 
$+W(I , J)*CH9) 

$-0. 8*(U( I , J)*CH2+W( I , J)*CH10)*SS/R(I)*S**2+4 .0*0M 
$*S**3*SS/R(I) 

AR(3 , 4)--4 ,O*OM*S-0 . 8 * ( U ( I , J ) *CH2+W( I , J ) *CH i 0 ) 

AR(3 , 6)--0 .8*(GW+CH8*U (I , J)+WCI , J)*CU9)-4 . 0*0M*C 
AR(3 , 7 )»-TD*2 .2/3.0 

AR(4,1)=--0.4*(U(I , J)*CH2+W(I , J)*CHlO)-2 .0*OM*S 
AR(4 ,3)— 0.4*(-SS/R(I)*S*C+SS/AL)*(GU+U(I » J)*CH5 
$+W(I , J)*CH4) 

$+0. 4*SS/R(I)*S**2*(GUS+U( I , J)*CH1+W(I , J)*CH3)+2 ,0 
$*OM*S*(l .0 

$+SS**2/R(I)**2*S**2-2 .0*SS**2/R(I )/ AL*S*C) 

AR(4 ,4)— 0.4*(GUS+U(I , J)*CH1+W(I, J)*CH3)-0.4*(-8S/ 
$R(I)*S*C 

$+SS/AL>+0 .4*SS/R(I)*3**2*(U(I , J)*CH2+W(I , J)*CHiU) 
$-2 .0*OM*C*(-SS/R(I)*S*C+SS/AL)-1 . 5*TD/TE 
AR(4 , 5)— 0 .4*(GW+U( I , J)*CH8+W( I . J)*CH9)-2 .0*OM*C 
AR(4,6)— 0.4*(GU+U(I , J)*CH5+W(I , J)*CH4)+2 .0*0M*S 
$*(-SS/R(I)*S*C 
$+SS/AL) 

AR(4 , 7 )-2 .2/3.0*TD*3S/R(I )*S**2 
AR(5, 1 )«-2 .0*0M*C*SS/R(I)*S**2 

AR ( 5 , 2 ) — 0 . 4 * ( GU+U ( 1 , J ) *CH5+W ( I , J ) *CH4 ) + 2 . 0 *OM* S * 
$(-SS/R(I)*S*C 
$+SS/AL) 

AR(5 ,4)-2 ,0*0M*C*( 1 .0+SS**2/R(I)**2*3**2 ) 

AR(5 , 5)»-0.4* (GUS+U(I , J)*CH1+W( I , J)*CH3 )-2 .0*OM 
$*S**3*SS/R(I) 

$+2 .0*OM*C*(-SS/R(I)*S*C+SS/AL)-1 .5*TD/TE 
AR(5 ,6)— 0.4*(-SS/R(I)*S*C+33/AL)*(GU+U(I , J)*CH5 
$+W(I,J)*CH4) 

$+2 .0*0M*S*(1 .0+SS**2/R(I)**2*S**2-2 .0*SS**2/R(I)/ 
$AL*S*C) 

AR(5 , 7)— TD*2 .2/3.0*(SS**2/R(I)**2*C*3**3-S8**2/ 
$R(I)/AL*S**2) 

AR(6,2)-“ 0.4*CGW+U(I, J)*CH8+W(I , J)*CH9)-2 .0*OM*C 
AR(6 , 3)-2 .0*0M*C*(1 . 0+SS**2 /R( I ) **2 *S **2 
$-2.0*SS**2/R(I)/AL*S*C) 

AR(6,4)— 2.0*OM*C*SS/R(I)*S**2 
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AR(6,5)— 0.4*(U(I . J)*CH2*rW(I,J)*CH10)-2.0*0M«;i 
AR(6,6) — 0 4 4*(-SS/R(I)*S*C+SS/AL)<»(GW+U(I, J)*CHb 
$+W(I,J)*CH9) 

$+0.4*SS/R(I)*S»»*2*(0(I, J)*CH2+W(I , J)*CH10)-2 .0>^0l-l 
$*C*(-SS/R(I) 

$*S*C+SS/AL)+2 4 0*OM*S*SS/R(I)*S**2+2.0*OM*C*(-JiS/R 
$(I)*S*C+;iS/AL) 

$-1 .5*TD/TE 

AR(6,7)-2.2/3.0*TD»SS/R(I)>»3**2 

DO 4930 ll-l ,3 

IJ-II+Jl 

AR(XI, ll)-i .0/AR(ll, 11) 

AR(I1, 7)-AR(Il, 7)*kAR(lI,ll) 

DO 493i N-IJ.6 
AR(11,N)»AR(1I,N)*AR(11, II) 

4931 CONTINUE 

DO 4932 N>1J,6 

AR(N,7)-AR(N,7)-AR(II,7)*AR(N.11) 

4932 CONTINUE 

DO 4933 N-IJ,6 
DO 4934 M-IJ,6 

AR(14,N)-AR(M,N)-AR<H»11)*AR(11.N) 

4934 CONTINUE 

4933 CONTINUE 
4930 CONTINUE 

AR(G, 7)-AR(6,7)/AR(6»6) 

NV«6 

4943 NV-NV-1 
IF(NV.EQ.O) GO TO 4944 
NX-NV+i 

DO 4930 N>NX,6 

AR(NV,7)-ARCNV,7)-ARCNV,N)»AR(N,7) 

4930 CONTINUE 

GO TO 4943 

4944 RS(1 ,J)-AR(1 ,7) 

RS(2,J)-AR(2,7) 

RS<3,J)-AR(3.7) 

RS(4, J)-ArU , 7) 

RS(5»J)-AR(5,7) 

RS(6, J)-AR(6,7) 

PRE— (GU*R3(5. J)+GW*RS(0,J) )-TD 

WRITE (6, 6330) I » J , GU , GW , RS ( 1 , J ) . RS( 2 « J ) , RS( 3 , J ) , RS 
$(4,J),RS(5,J),RS(6,J),PRE 
4905 CONTINUE 

DO 8921 N9-1 ,6 
RS(N9,1)«RS(N9»2) 

RS(N9, 1MM)>RS(N9,2) 

8921 CONTINUE 
MM-0 

4013 MM-MH+1 
NEW S 

DO 4020 J«1,IMM 
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RUV-0.0, 

DN-DNO/4.0, 

IF(J.LT.IMMi.0R.J.GT.IMM2) DN-DNO , 

WCOV»-SS/R(I),*:i**2*UCl , J) + (-SS/R(I)*S*C+SS/AL)*V(I ,J) 
$+( 1 .0+SS**2/R(I)**2*S**2-2.0*SJ**2/R(i)/AL*S*C)*W(I, J) 
UCOV-U(I , J)-S3/RU)*S**2*WCI, J) 

VCOV-V(I , J)+(-SS/R(I)*S*C+SS/AL)*W(I , J) 

IF( CJ“ID) .LT . 1) GO TO A145 
IF( (J-ID) .GT.IMML) GO TO 4146 
UNN=U(1-1 , J-ID) 

FNN“F(l-i , J-ID) 

GO TO 4141 

4145 UNN«U(I-1,1) 

PNN»P(1-1,1) 

GO TO 4141 

4146 UNN=U(I-1 , IMML) 

PNN=P(I-1 ,IMHL) 

4141 DN=DN0 

IF( J .GE . IMMl .AND. J. LE. IMM2)DN-DN0/4.0 

A(i , J)=-V(l, J)/2.0/DN 

A(2 , J)=‘(UN(I , J)-UL( I , j; )/2 .O/OS 

AC3,J)=VCi,J)/2.0/UN 

IF(J.EQ. 1 .OR. J.EQ. HIM) GO TO 9011 

PSN=PCl,J+l)-P(i,J-l> 

RSN = RS(5 , J+1 )-RS(5 , J-1 ) 

GO TO 9012 

9011 RSN*0.0 
PSN=0 .0 

9012 A(4, j;=-U(I, J)**2-W(I ,J)*(U(I+1 , J )-UNN >/2 .0/DR 
$-W(I, J)**2^CH1 I 

$-2.0*U(l,J)*W( I , J)>*'CH3-2.0*U(1 , J ) «V ( I , J ) *0H5-2 . 0*V 
$(I , J >*W( I , J)*CH4 
$ + 2 .0*OM’^S*WC0V-RSN/2 .0/DN 

$+1.0/DENS*C(l .0+SS**2/R(I)**2*S**4 )*(FN(I,J) 

$-PL(I , J) )/2 .0/DS 

$+(SS*2/R(I)**2*S**3*C-SS**2/R(I)/AE*S**2 )*PSN 
$/2 . 0/DN+SS/R(I)*S**2*(P(I+l , J ) -PNN ) / 2 . 0 /DR ) 

IF( J . EQ. IMMl ) GO TO 4019 
IF(J.EQ. IMM2) GO TO 4018 
GO TO 4020 

4019 A( 1 , J)=-V(I ,J)/8.0/DN 

A(2 , J)*(UN(I , J)-UL(I , J) )/2.0/DS-3.0/8.0*V(I ,J)/DN 

A(3,J)*V(I,J)/2.0/DN 

GO TO 4020 

4018 A( 1 , J)=-V(I , J)/2.0/DN 
A(3 ,J)=V(I,J)/8.0/DN 

A(2 , J)=3 .0/8.0*V(I, J)/DN+(UN(I , J)-UL(1 ,J))/2.0/DS 

4020 CONTINUE 

A(4, 1)»A(4,1)-A(1 ,1)*U(I,1) 

A(4 , IMM) = A(4 »IMM)-A(3, IMM)*U(1 , IMM) 

CALL ELIM(A,IMM) 

DO 4101 J-1 , IMM 

U(I, J)-U(I, J)+SOR*(A<4, J)-U(I, J)> 

4101 CONTINUE 


n a 


158 


NEW R 
C 

DO 4030 
RWV-0.0 
ON-DNO/4.0 

IF(J.LT.lMMi .OR. J.GT.IMM2) DN«DN0 

WCOV--SS/R(I)*S**2*U(I » J>+(-SS/R(l)*S»'C+SS/AL)*V(l , J) 
$+( 1 .0+3S**2/R(I)**2*S**2-2 .0*SS**2/R(l)/AL*a*C)*W( I , J) 
UCOV-U(I , J)-Sb/R(l)*S**2*W(I, J) 

VCOV-V(l, J)+(-SS/R(I)*3*C+SS/AL)*W(l, J) 

IF( (J-ID) .LT. 1 ) GO TO 4165 
IF( (J-ID) .GT . IMML)GO TO 4166 
WNN» W(I-1 , J-ID) 

PNN-F(I-1 , J-ID) 

GO TO 4161 

4165 WNN- W(I-1 ,1) 

PNN=P(I-1 , 1 ) 

GO TO 4161 

4166 WNN- W(l-1 , Il'lML) 

PNN-P(I-1 , IMML) 

4161 A(1 , J) — V(l, J)/2.0/DN 

A(2,J)-(W(I+1 , J)-WNN)/2.0/DR 
A(3,J)-V(1,J)/2.0/DN 
IF(J.EQ.1 .OR.J.EQ.IMM) GO TO 9111 
RSN-RS(b , J+1 )-RS(6 , J-1 ) 

PSN»P(I,J+1)-P(I,J-1) 

GO TO 9112 

9111 PSN-0.0 
R3N-0.0 

9112 A(4, J)— U(I , J)*(WN(I , J)-WL(I , J) )/2 .O/DS-2 .0*OM*C*VCOV 
§-2 .0*OM*UGOV*S-W(I» J)*’^2*Clii2-U(2, J)**2»'CH2-2 .0* 

$U(I, J)*W(I ,J)*CH10 

$-2 ,0*W(I , J)*V(I-, J)*CH9-2 .0*0(1 , J)*V(I , J)*CHti 
$-R3N/2 .0/DN-RS( 1 , J ) *CH2 -RS ( 5 , J ) *CH8 
$-1.0/DEHS*(SS/R(i)*S**2*(PN(l,J)-PL(I*J))/2.0/DS 
$+(SS/R(I) 

$*S*C-aS/AL)*PSN/2 .0/DN+(P(I+l , J)-PNN) 

$/2.0/DR) 

IF(J.EQ.IMMl) GO TO 4029 
IF(J.EQ. IMM2) GO TO 4028 
GO TO 4030 

4029 A(1 , J) — V(I,J)/8.0/DN 
A(3 , J)-V(I , J)/2.0/DN 

A(2 , J)-(W(I+1 , J)-WNN)/RD/2 .0-3 .0/8.0*V(I, J)/DN 
GO TO 4030 

4028 A( 1 , J) — V(I , J)/2 .0/DN 
A(3,J)-V(I ,J)/8.0/DN 

A<2 , J)-3 .0/8.0*V(I , J)/DN+( W( 1+1 , J ) - WNN ) / 2 . 0/DR 

4030 CONTINUE 

A(4,l)«A(4.1)-A(l,l)*W(i,l) 

A(4 , IMM)-A(4 , 1MM)-A(3, IMM)*W(I , IMM) 

CALL ELIM(A,IMM) 

DO 4102 J-1, IMM 
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W(I,J)-W(I,J)+S0R*(A(4, J)-W(I,J)) 

4102 CONTINUE 
UMl-200.0 

DO 4103 J-1 , IMM 

IFCJ.LT.IMMl .OR. J.GT.IMM2)GO TO 4103 
IF(Ud.J) .GE.UMDGO TO 4103 
UM1-U(I,J) 

NT-J 

4103 CONTINUE 
NEW N 


DO 60 30 j-1/xmi 
RWV“0 .0 
DN«DNO/4.0 

If (J . LT. IMMl .OR. J.GT. IMM2) DN^DNO 

WCOV=*-SS/R(I)*0^*2*U(I , J) + (^yS/R(I)*3*C+SS/AL)‘*'V(I , J) 

$ + ( 1 .0+3S**2/R(I)**2*3**2-2 .0*S£»**2/R(I)/AL*S*C)*W(I , J) 
UC0V“U(I,J)-3S/R(I)*S**2*W(I, J) 

VCOV=V(I ,J)+(-SS/R(I)*S*C+SS/AL)*W( I . J) 

IF( (J-ID) .LT. 1 ) GO TO 6165 
IF( (J-ID) .GT . IMMDGO TO 6166 
VNN= V(I-1,J“ID) 

FNN»P( I-l , J-ID) 

GO TO 6161 

6165 VNN= V(l-1 , 1 ) 

PNN = P( I-l , 1 ) 

GO TO 6161 

6166 VNN= V(I~1 , IMML) 

PNN=P(I-1 , IMML) 

6161 AC i , J)“-V (I , .0/DN 
A(2,J)“U(I,J)/0S 
A(3.J)*V(1,J)/2.0/DN 
IFCJ.EQ. 1 .OR. J.EQ. IMM) GO TO 9211 
PSN=P(I , J+1 )-P(I , J-1 ) 

GO TO 9212 

9211 PSN=O.U 

92 12 A(4 , J) = U( I , J)*VLL(I , J)/IJS-W(I , J)*(V(1-H , J ) - VNN ) / 2 . 0 / DR 
$-V(I , J)**2*CH6-W(1 , J)**2*CH14-2.0*M(I, J)*CH15 
$-2 .0*V(I , J) *W(I, J>*CH16+2 .0*OM*WCOV*C 

$-1 .O/DENS’^C (SS*»''2/R(I)**2*S**3*C-SS**2/R(I)/AL*S**2 ) 
$*(PN(I, , J))/2.0/DS + ( 1 .0+SS**2/R(1)**2*(S*C)**2 

$-2 .0*SS**2/R(I)/AL'^b’<*C)*PSN/2.0/DN 
$+(SS/R(I )*S*C-SS/AL)*(P(I + 1 , J)“I*NN)/2 .0/DR) 

IFCJ.EQ. IMMl ) GO TO 6029 
IFCJ.EQ.IMM2) go TO 6026 
GO TO 6030 

6029 ACl »J)“-V(I,J)/6.0/DN 
A(3 , J)“VC 1 , J)/2 .O/DN 
AC2»J)-U(I»J)/DS-3.0/3.0*V(I,J)/DN 
GO TO 6030 

6026 ACl , J)-~V(I,J)/2.0/DN 
AC3,J)*VCI , J)/8.0/DN 
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A(2 * J)«3.0/i).0*V(I , J)/DN+U(I , J)/DS 
6030 CONTINUi: 

A(4,1)»A(4,1)-A(1»1)*V(I,1) 

A(4 , IMM)=A(4 , IMM)~A(3, IMM)*V(I , IMM) 

CALL IiLIM(A,IMM) 

DO 6102 J=-i,IMM 

V(l, J)»V(I,J)+SOR*(A(4,J)~V(I, j;) 

6102 CONTINUt 

NEW P 

DO 4320 
RUV=0 . 0 
DN=DN0/4 .0 

IF ( J . LT . IMM 1 . OR . J . GT . IMM2 ) DN=DN0 

WCOV=-SS/R(i)*S**2*U(I , J) + (-ilS/R(I)*S*^C+SS/AL)*V (i , J) 

$ + (l.U+SS*’=f2/R(I)**2*S**2-2.0«CS**2/R(I)/AL*G’^C)*W(I,J) 
UCOV=U(I , J)-SE/R(I)*S**2*W(I, J) 

VC0V=V(I, J)+(-3S/R(l)*S«C+SS/AL)*W(I, J) 

IF( (J-ID) .LT. 1 ) GO TO 4545 
IF( (J-ID) .GT. IMiML) GO TO 4546 
UNN=U(1-1 , J-ID j 
PNN=P(1-1 , J-ID) 

VNN=V(I-i , J-ID) 

WNN=W(I-1 , J-ID) 

UNNL=UL(I-1 , J-ID) 

UNNN=UN(I-1 , J-ID) 

VNNL=VLL(I-1 , J-ID) 

VNNN=VN(I-i , J-ID) 

WNNL=WL(I-1 , J-ID) 

WNNN=WN(I-1 , J-ID) 

PNNN=PN(I-1 , J-ID) 

PNNL=PL(I-1 , J-ID) 

GO TO 4541 
4543 0NN = a(l-l , 1 ) 

PNN = P(I-1 , 1 ) 

VNN = VCI-1 . 1 ) 

WNN=W( I-l , 1 ) 

UNNN = UN(l-i , 1 ) 

UNNL = aL(I-l , 1 ) 

VNNN = VN(I-1 , 1 ) 

VNNL = VLL( 1-1 , 1 ) 

PNNN=PN(I-1 , 1 ) 

PNNL = PL(I-1 , 1 ) 

WNNN = WN(I-1 , 1 ) 

WNNL = WL(I-1 , 1 ) 

GO TO 4541 

4546 UNN=U( X-1 , IHML) 

PNN=P(I-1 , IMML) 

VNN=V(I-i , IMML) 

WNN=W( I-l , IMML) 

UNNN=UN( I-l , IMML) 

UNNL=UL(I-1 . IMML) 

VNNN=VN(I-1 , IMML) 
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VNNL“VLL(I-1 , IMML> 

FNNN-FNU-l , IMML) 

PNNL-PLCl-1 , IMliL) 
wWNN-WN(I“l , IMML) 

WNI” »WL(I-1 . IMML) 

4341 DN-bi'lC 

FSl-1 .0+SS**2/R(I)**2*S**4 

PS2=-SS**2/R(I)**2*S**3*C-SS**2/R(I)/AL*S**2 

PS3*SS/R(I)*S**2 

PS4»1 .0+SS**2/R(I)**2*(S*C)**2-2.0*SS**2/R(I)/AL*S*C 

PS3-SS/R(I)*S*C-SS/AL 

PSO = i .0 

IFCJ.GE . IMlli .AND. J .LE. Il’lM2li)N»DN0/4 .0 
A( L , J)“l .0/DN**2 
A( 2 , J)»-2 .0/DN‘^*2 
A(3 , J)-l .0/DN**2 

IF(J.EQ.l.OR.J.EQ.XMii) GO TO 9311 
PSN1=PN(I , J+i 1+PL(I , J“1 )-PNU , J-1 )-PL(I . J+1) 

PSN2-P(i+1 , J+1 )-P( I+i , J-i > 

RSN=RSC2 , J+i)-RSC2 , J-1) 

WSN = W(I , J+1)*W(1 , J-1 ) 

GO TO 9312 

9311 PSN1=0.0 
PEN2“0,0 
WSN“0.0 
RSN=0.0 

9312 A(4, J)— PS1*(PN(I , J)+PL(I »J)«2.0*P(1,J))/DS**2 
$-PS6*(P(I + l »J)+PNN“2.0*P(X ,J) )/DR**2-0.3*PS2* 
$PSN1/DN/DS“0 .5*PS3*(PN(I+1 , J) 

$+PNNL-PNNN“PU(l+l . J) )/DR/DS 

IF( (J-ID) .LT . 1 .OR. ( J-ID) .GT . IMML) GO TO 4707 
PNNl=P(I-l , J-i )-P(I-l , J+1 ) 

GO TO 4708 
4707 PNN1*0.0 

4 7 08 A(4 , J)=A(4 . J)-0 .5*PS5*(PSN2+PNN1 )/DR/DN 

$+DKNS*(( (UN(I,J)-UL(I.J))/2.0/DS+U(l , J)*CUi+W(I , 
$J)*CH3*C 

$(W(I+1 , J)-WNN)/2.0/DR+U(I,J)*CH3+W(I ,J)*GH11) 

$-( (U(I , J)-UNNl/2 .0/DR+U(I . J>*CU3+W( I , J)*CUii )*C(WN(1 » Jl 
$-WL(I,J))/2.0/DS+U(I,J)*CH2+W(I,J)*CH10)- 
$(RSN-2.0*RS(2, J) )/DN**2 

$-2.0*(-OM*S*((WN(I , J)-WL(I , J>)/2.0/DS“(U(I+l , J)-UNN 
$)/2 . 0/DR)+0M*C*( (V(I+1 . J)-VNN)/2.0/DR-W3N 
$/2.0/DN))) 

IFCJ.Eq.TMMl ) GO TO 4519 
IF(J.EQ. IMM2) GO TO 4518 
GO TO 4520 

4519 A( I , J)-0. 10/DN**2 
A(2, J) — 0.50/DN**2 
A(3 , J)-0.4/DN**2 
GO TO 4520 

4518 A(1 , J)-0.4/DN**2 

A(3 , J)-0. 100/DN**2 
A(2, J)“-0.50/DN**2 
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4520 CONTINUE 

A(4,i)=A(4,l)~A(l,l)*P(I,l) 

A(4 , IMM)=A(4 , XIiM)-A(3 , IMM) ( I , IMM) 

CALL ELIM(A,IMM) 

DO 4501 

P(1,J)=P(I, J)+SOR*(A(4, J)-P(I*J)) 

4501 CONTINUE 

IF(MM.GE.2) GO TO 4600 
GO TO 4015 

4600 DO 41020 J«1 , TOI 

IF((J-ID) .LT.l) GO TO 41145 
IF( (J-ID) .GT, IHML) GO TO 41146 
ENNN= EN( 1-1 , J-ID) 

DINN= DI(I-1 , J-ID) 

GO TO 41141 

41145 ENNN = EN(I-1 , 1 ) 

DINN= DI( I-l , 1) 

GO TO 41141 

41146 ENNN=£N(I-1 , IMHL) 

DINN^ DI(I-1 , IMIIL) 

41141 DN=DN0 

If (J.GE.IMMi .AND.J.LE.IHM2) DN=DN0/4.0 

Q1=0.0 

Q2=0 .0 

Q3=0.U 

Q4 = 0 . 0 

If (J,£Q. 1 .OR. J.EQ. IMM) GO TO 4077 
IF( J .EQ. IMMl ) GO TO 4075 
IF(J.EQ. IHM2) GO TO 4076 

Q1=U(I , J) *CH5i-W(I , J)*CH4 + (U(I , J+i )-U(I , J-1 ) )/2 .0/DN 
Q2=U(I ,J)*CHd+W(I , J)*CH9+(W(I , J+1 ) -W( I . J-1 ) ) /2 . 0/DN 
Q3=(UN(I ,J)-UL(1 , J) )/2.0/DS+U(I , J)*CHH-W(I,J)*CH3 
Q4=(WN(I , J)-UL(I , J) )/2.0/DS+U(I , J)*CU2+W(I , J)*CH10 
GO TO 4077 

4075 Q1=U.5/DN*(U(I , J+1 )-0. 75*U(I , J)-0.25*U(I , J-1 )) 

$+U(I , J)*CH5+W(I , J)*CU9 

02=O.5/DN*(W(I , J+i )-0. 75*W(I , J)-0.25*W(I , J-1 )) 

$-t-U( I ,J)*Ciia+W( I , J)*CH9 
GO TO 4077 

4076 Q1=0. 5/DN*(0.25*U(I , J+i )+0. 75*U(I, J)-U(I , J-1 ) ) 

$+U(I , J)*CH5-f-W(I , J)*CH4 

■Q2=0 . 5/DN« (0 . 25*W(I , J+1 )+0 . 75*W(I , J)-W(I , J-1 ) ) 

$+U(I , J)*CH8+W( I , J)*CH9 

4077 E=0 . 1*RS(2 , J )*EN( I , J)/DI(I , J) 

PR=-RS(1 , J)*Q3-RS (5 ,J)»'Q 1-RS (4, J)*Q4 -RS (6, J)*Q2 

IF( J . LQ. IMl'II ) GO TO 41019 

IF( J . EQ . IMM2 ) GO TO 410lo 

A(ijJ)=-VCl,J)/2.0/DN-E/DN*’^'2 

A(2,J) = U(I,J)/DS+W(I,J)/DR +£*^2. 0/UN* *2 

A(3,J)=V(I,J)/2.0/DN-E/DN**2 

A(4 , J)=PR-DI( i , J)+U(I , J)*EN1 (J)/DS+W(I , J)*ENNN/DR 
G( 1 ,J) = -V(I , J)/2 .0/DN-E/l .3/DN"*2 
G(2 ,J)=U( I , J)/DS+W( I , J)/DR+E*2 .0/1 .3/DN**2 
G(3 , J) = V(I , J)/2 .0/DN-E/ 1 . 3/DN**2 
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G(4, J)=U(I,J)*DI1(1)/DS+ W( I,J)*DINN/DR 
$+i .44*DI(I, J)/EN( I , J)*PR 
$-1 .92*DI(I,J)**2/EN(I , J) 

GO TO 41020 

41019 A(1 , J)— V(I, J)/8.0/UN- E/10.0/DN**2 

A(2,J)»U(I,J)/DS-3.0/8.0*V ( I , J ) /DN+WN ( I , J ) /DR 

$+0,5*E/DN**2 

A(3 , J)=V(I, J)/2.0/DN-E/2.5/DN**2 

A(4 , J)-PR-DI (I » J)+U(I , J)*EN1 (J)/DS+W ( I , J ) *ENNN / DR 
G(1 , J)«-V(I, J)/8 .0/DN-E/13.0/DN**2 
G(2,J)»U(I, J)/DS~3.0/8.0*V(I, J)/DN 
$+W(I , J)/DR+1 .0/2 .6*E/DN**2 
G(3 , J)=V(I, J)/2.0/DN-'E/2.5/l .3/DN**2 
G(4 , J)=U(I, J)*DI1 (J)/DS+W (I , J)*DINN/DR 
$+1 .44*DI(I , J) /EN(I , J)*PR 
$ -1 . 92*DI( I , J)**2/EN( I , J) 

GO TO 41020 

41018 A(1 , J)=-V(I , J)/2.0/DN-E/2 ,5/DN**2 
A(3 , J)»V( I , J)/8 .0/DN-E/10.0/DN**2 
A(2 ,J)=U(I,J)/DS+3.0/8.0*V ( I , J ) / DN+WN (I , J ) / DR 

$+0.5*E/DN**2 

A(4,J)*U(I , J)*EN1 C J)/DS+W(I , J)*ENNN/DR+PR-DI(I , J) 
G(1 ,J)=-V(I,J)/2.0/DN-E/2.5/l .3/DN**2 
G(3 , J)=V(I , J)/8 .0/DN-E/l 3.0/DN**2 
G(2, J)=U(I ,J)/DS+3 .0/8.0*V(I, J)/DN+W(I , J)/DR 
$+0.5/1 .3*E/DN**2 

G(4 , J)=U(I » J)*DI1 (J)/DS+W(I , J)*DINN/DR 
$+1 .44*DI(I , J)/EN(I , J)*PR 
$-1.92*DI(I,J)**2/EN(I.J) 

41020 CONTINUE 

A(4, l)«A(4, 1)-A(1 . 1)*EN(I, 1) 

A(4 , IMM)=A(4 , IMM)-A(3 , IMM)*EN( I , IMM) 
G(4,l)=G(4,l)-G(l,l)*DI(I.l) 

G(4 , IMM)=G(4 , IMM)-G(3 , IKM)*DI(I 
CALL ELIM(A,1MM) 

DO 41101 J=1,IMM 

EN(I ,J)=EN(I , J)+SOR*(A(4 , J)-EN(I , J) ) 

41101 CONTINUE 

CALL ELIM(G,IMM) 

DO 41107 J=1 , IMM 

Did , J) = DI(I , J) + S0R*(G(4, J)-DI(I, J) ) 

41107 CONTINUE 

C 

C 

Q itith ititit hitlt A 4c A 

C 

c 

N20-N101 (I) 

WRITE(6 , 2370)K, I 

WRITE (6, 2^70) (N10(I,J).U(I,N10(I,J)).V(1,N10 

$(I.J))»W(I,N10 

X(I,J)),EN(I,N10(I,J) ) ,DI(I,N10(I»J)) . .N20) 

4010 CONTINUE 

DO 4820 I»l,7 
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KiI«IMAX(I) 

WiUTii(i4,i>30) (U(I.J;.\!(I,J),W(I,J>,F(I»j;,iiN(X,J 
$),D1(I.J) ,J=i,KM) 

UL(I ,J)^UC1,J) 

VLLCI » J)**V(I , J) 

PL(I.J)«i»(I,J) 

£NL(I ,J>-KN(I , J) 

DILU.J)-DI(X ,J) 

4820 CONTINUE 
4700 CONTlNUt 
44444 STOP 
£NU 

SUUROUTlHt: ELIM(S.NN) 

DIMENSION S’( 4,300) 

NNOl-NN-1 

DO 1 11=1, NNOl 

S(2 , ID-1 .0/S(2, II) 

S(3 , II)-S(3 , XI)*S(2 , II) 

S(4,II)-SC4,II)*S(2,II) 

S(4, II+D-S(4,II+1)-S(4, 1I)*S(1 ,11+1 ) 

S(2, II+l )»S(2 , II+1)'S(1 , II+1)*S(3,II) 

1 CONTINUE 

S(4 ,NN)»S(4 ,NN)/S(2 ,NN) 

NV-NN 

3 NV-NV-i 

IF(NV.EO*0) GO TO 4 
NVAl-NV+1 

S(4 .NV)»S(4,NV)-S(3,NV)*S(4 ,NV + 1 ) 

GO TO 3 

4 RETURN 
END 
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